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INTRODUCTION 


This final report summarizes the research accomplishments under Contract NAG 
2-648 between NASA Ames Research Center and the Pennsylvania State University. 
More specific information on each topic covered in this report is contained in the 
manuals and preprints which have been submitted along with this report, as described in 
the Deliverables section of the Statement of Work for the subject contract. Also 
submitted with this report are 4 different FDTD computer codes and companion RCS 
conversion codes on magnetic media. In the remainder of this report the preprints, the 
computer codes and their user’s manuals will be summarized. 


COMPUTER CODES DELIVERED 

Under this effort a single three dimensional dispersive FDTD code for both 
dispersive dielectric and magnetic materials was to be developed and delivered, along 
with a User’s Manual. This code is included with this report on Magnetic Media, and is 
named "FDTDD", standing for the "D" version of a set of three-dimensional FDTD codes 
developed at Penn State. This code has the capability to calculate electromagnetic field 
interactions with objects which include both dispersive and non-dispersive dielectric and 
magnetic materials. Using a companion code, "RCS3D", the output from FDTDD can 
be converted to radar cross section vs frequency. A User’s Manual describing the theory 
and use of the computer code FDTDD and showing validation results is included with 
this final report as attachment [1]. 

In addition to FDTDD, simpler (and somewhat faster) three-dimensional FDTD 
codes which have more limited or no dispersive material capability, "FDTDA" through 
"FDTDC", are also being delivered along with this final report. Version "A" is for 
frequency-independent materials, version "B" for frequency-dependent dielectric (non- 
magnetic) materials, and version "C" for frequency-independent dielectric and magnetic 
materials. While version "D" includes all of these cases, the simpler versions "A" through 
"C" are somewhat easier to use and will run faster than the more complicated "D" 
version. They are therefore preferred when the more extensive capabilities of the "D" 
version are not needed. All versions use the same "RCS3D" companion computer code 
to convert time domain output to RCS vs frequency. The User’s Manuals for all four 
versions of the Penn State FDTD code are included with this report as attachments [1-4]. 


While three-dimensional FDTD codes are generally more useful, two-dimensional 
codes can be applied to geometries of interest, and have the advantage of requiring 
significantly less computer resources. A pair of two-dimensional FDTD codes, TEA" 
and "TMA", for transverse electric and transverse magnetic excitation, respectively, are 
also included with this final report on magnetic media. Also included are companion 
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computer codes "SWTEA" and "SWTMA", which convert the time domain output of 
"TEA" and "TMA" to scattering width vs frequency. The user’s manual for these codes is 
also included with this report as an attachment [5]. A theoretical description and 
validation results for the two-dimensional computer codes are contained in a preprint 
included with this report as attachment [6]. 


EXTENSIONS TO DISPERSIVE FDTD CAPABILITIES 

During this effort two extensions to our dispersive material FDTD capability were 
made. One was the extension to include magnetic materials which have more 
complicated frequency dependence of their permittivity than considered previously. The 
time domain magnetic susceptibility has the form of a damped sinusoidal wave, resulting 
from the oscillations of the microscopic currents in the material. Dispersive FDTD was 
extended to include this class of materials, and these results are reported in a PhD 
dissertation by Forrest Hunsberger that is included with this final report as attachment 
[7]. Some of this material has already been presented at the recent IEEE AP-S 
conference [8], and Dr. Hunsberger is currently writing papers based on this dissertation 
for journal submission. 

The other extension is to a time domain surface "impedance" formulation for lossy 
conductors. As explained in attachment [9], this removes a difficulty in applying FDTD 
to targets containing lossy materials. The difficulty is due to the necessity of reducing 
the size of the FDTD cells which are inside the conducting materials. With the FDTD 
surface "impedance" described in [9], FDTD cells need not be located inside the 
conducting material but only at the surface. The attachment [9] preprint has already 
been accepted for publication in IEEE Trans, on Antennas and Propagation. These 
results have in part been presented at the recent IEEE AP-S symposium [10]. 

The approach used in [9] also extends the basic frequency dependent FDTD to a 
wider class of materials. This is due to the method used in [9] to obtain the necessary 
time domain convolution coefficients. With this method ANY material whose frequency 
dependent behavior is known can have this information be transformed to the time 
domain, and then with application of Prony’s method (see [9]) the necessary exponential 
coefficients for applying dispersive FDTD can be obtained. 


PLATE SCATTERING 

During the course of this research effort three different plate scattering 
geometries were given special consideration. The first of these was suggested by Dr. 
Randy Jost. It involved a conducting plate coated with a dielectric layer. This geometry 
cannot readily be modeled using other scattering calculation methods. Scattering results 
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for this geometry were calculated using one of our FDTD codes (the "FDTDA" code 
furnished under this effort). The results were documented and are attached to this report 
[11]. The results were previously furnished to Dr. Jost, and according to Dr. Jost they 
showed significant correlation with measurements. We hope to continue this 
investigation. 

The second plate geometry considered was the "business card" plate geometry 
being used as a test case by Dr. Alex Woo of NASA Ames. This geometry is challenging 
since edge waves contribute to the scattering. Moment Method codes require 
approximately 10 times the usual density of modes to accurately compute the scattering 
from this geometry. Our results, shown in attachment [12] of this report, shown the 
difference in RCS computed using FDTD and a Moment Method code. Dr. Woo 
indicated that our results were quite accurate, considering that we used only 10 cells per 
wavelength running on a 486 PC. We hope to run this data again using more cells per 
wavelength, and compare our results directly with measurements. 

The third topic considered was extending FDTD to include modeling of thin 
impedance sheets. While applicable to arbitrary shapes and in both two and three 
dimensions, the extension was validated by calculating RCS from thin flat impedance 
sheet plates. The results are reported in attachment [13]. 


NONLINEAR MATERIALS 

Having extended the FDTD method to dispersive materials, another class of 
materials which it would be desirable to include are nonlinear materials. During this 
effort the current flowing in a wire antenna loaded with a nonlinear diode was 
successfully computed using FDTD. The results are shown in the attached preprint [14]. 
With this capability RCS from scatters including nonlinear loads or bulk materials may 
be computed using FDTD. Further extensions could include scattering of short pulses 
from nonlinear dispersive materials including ferrite absorbers. 


WIRE ANTENNAS 

The accuracy available from the FDTD method was demonstrated by computing 
the self impedance of a wire dipole antenna and the mutual impedance between two 
wire antennas. These results are reported in attachment [15]. 
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CONCLUSIONS 


During this effort the tasks specified in the Statement of Work have been 
successfully completed. The extension of FDTD to more complicated materials has been 
made. A three-dimensional FDTD code capable of modeling interactions with both 
dispersive dielectric and magnetic materials has been written, validated, and documented. 
This code is efficient and is capable of modeling interesting targets using a modest 
computer work station platform. 

However, in addition to the tasks in the Statement of Work, a significant number 
of other FDTD extensions and calculations have been made. RCS results for two 
different plate geometries have been reported. The FDTD method has been extended to 
computing far zone time domain results in two dimensions. Finally, the capability to 
model nonlinear materials has been incorporated into FDTD and validated. 

The FDTD computer codes developed have been supplied, along with 
documentation, and preprints describing the other FDTD advances have been included 
with this report as attachments. 
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I . INTRODUCTION 

The Penn State Finite Difference Time Domain Electromagnetic 
Scattering Code Version A is a three dimensional numerical 
electromagnetic scattering code based upon the Finite Difference 
Time Domain Technigue (FDTD) . The supplied version of the code 
is one version of our current three dimensional FDTD code set. 
This manual provides a description of the code and corresponding 
results for the default scattering problem. The manual is 
organized into fourteen sections: introduction, description of 

the FDTD method, Operation, resource requirements, Version A code 
capabilities, a brief description of the default scattering 
geometry, a brief description of each subroutine, a description 
of the include file (COMMONA. FOR) , a section briefly discussing 
Radar Cross Section (RCS) computations, a section discussing the 
scattering results, a sample problem setup section, a new problem 
checklist, references and figure titles. 

II. FDTD METHOD 

The Finite Difference Time Domain (FDTD) technique models 
transient electromagnetic scattering and interactions with 
objects of arbitrary shape and/or material composition. The 
technique was first proposed by Yee [ 1 ] for isotropic, non- 
dispersive materials in 1966 ; and has matured within the past 
twenty years into a robust and efficient computational method. 

The present FDTD technique is capabable of transient 
electromagnetic interactions with objects of arbitrary and 
complicated geometrical shape and material composition over a 
large band of frequencies. This technique has recently been 
extended to include dispersive dielectric materials, chiral 
materials and plasmas. 

In the FDTD method. Maxwell's curl equations are discretized 
in time and space and all derivatives (temporal and spatial) are 
approximated by central differences. The electric and magnetic 
fields are interleaved in space and time and are updated in a 
second-order accurate leapfrog scheme. The computational space 
is divided into cells with the electric fields located on the 
edges and the magnetic fields on the faces (see Figure 1 ) . FDTD 
objects are defined by specifying dielectric and/or magnetic 
material parameters at electric and/or magnetic field locations. 

Two basic implementations of the FDTD method are widely used 
for electromagnetic analysis: total field formalism and scattered 
field formalism. In the total field formalism, the electric and 
magnetic field are updated based upon the material type present 
at each spatial location. In the scattered field formalism, the 
incident waveform is defined analytically and the scattered field 
is coupled to the incident field through the different material 
types. For the incident field, any waveform, angle of incidence 
and polarization is possible. The separation of the incident and 
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scattered fields conveniently allows an absorbing boundary to be 
employed at the extremities of the discretized problem space to 
absorb the scattered fields. 

This code is a scattered field code, and the total E and H 
fields may be found by combining the incident and scattered 
fields. Any type of field quantity (incident, scattered, or 
total) , Poynting vector or current are available anywhere within 
the computational space. These fields, incident, scattered and 
total, may be found within, on or about the interaction object 
placed in the problem space. By using a near to far field 
transformation, far fields can be determined from the near fields 
within the problem space thereby affording radiation patterns and 
RCS values. The accuracy of these calculations is typically 
within a dB of analytic solutions for dielectric and magnetic 
sphere scattering. Further improvements are expected as better 
absorbing boundary conditions are developed and incorporated. 

III. OPERATION 


Typically, a truncated Gaussian incident waveform is used to 
excite the system being modeled, however certain code versions 
also provide a smooth cosine waveform for convenience in modeling 
dispersive materials. The interaction object is defined in the 
discretized problem space with arrays at each cell location 
created by the discretization. All three dielectric material 
types for E field components within a cell can be individually 
specified by the arrays IDONE(I, J,K) , IDTWO(I, J,K) , 

IDTHRE(I, J,K) . This models arbitrary dielectric materials with n 
= Mo* By a n obvious extension to six arrays, magnetic materials 
with n f can be modeled. 


Scattering occurs when the incident wave, marched forward in 
time in small steps set by the Courant stability condition, 
reaches the interaction object. Here a scattered wave must 
appear along with the incident wave so that the Maxwell equations 
are satisfied. If the material is a perfectly conductive metal 
then only the well known boundary condition 


E 


scat 

tan 


inc 

~Btan 


(1) 


must be satisfied. For a nondispersive dielectric the 
requirement is that the total field must satisfy the Maxwell 
equations in the material: 


Vx E tot 


Vx ( E inc + E scat 


1 3H tot 
M 0 dt 


1 d (H inc +H scat ) 

M 0 


dt 


(2) 
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VxH = Vx (H + H ) =e 


8 E 


tot 


rtot n„ /trine , TT scat^ MjE COt 


at 


(3) 


=e a(E inc + E sc ^M +g(E inc +E scat } 


at 


(4) 


Additionally the incident wave, defined as moving unimpeded 
through a vacuum in the problem space, satisfies everywhere in 
the problem the Maxwell equations for free space 


VxE inc 


i aH inc 
at 


(5) 


VxH inc 


= e 


aE inc 

at 


( 6 ) 


Subtracting the second set of equations from the first yields the 
Maxwell equations governing the scattered fields in the material: 


Vx E scat 


i aH scat 
M 0 dt 


(7) 


aE 


inc 


Vxh = (e -e ) 

0 at 


+aE ,nc + e 


dE : 


scat 


at 


+ (XE 


scat 


(8) 


Outside the material this simplifies to: 


Vx E 


scat 


1 dH 


scat 


M 0 


(9) 


VxH scat = e 


aE 


scat 


at 


(10) 


Magnetic materials, dispersive effects, non-linearities, 
etc., are further generalizations of the above approach. Based 
on the value of the material type, the subroutines for 
calculating scattered E and H field components branch to the 
appropriate expression for that scattered field component and 
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that component is advanced in time according to the selected 
algorithm. As many materials can be modeled as desired, the 
number equals the dimension selected for the flags. If materials 
with behavior different from those described above must be 
modeled, then after the appropriate algorithm is found, the 
code's branching structure allows easy incorporation of the new 
behavior. 

IV. RESOURCE REQUIREMENTS 

The number of cells the problem space is divided into times 
the six components per cell set the problem space storage 
requirements 

Storage = NC x 6 components/cell x 4 bytes /component (H) 

and the computational cost 

Operations = NC x 6 comp/cell x 10 ops/component x N (12) 

where N is the number of time steps desired. 

N typically is on the order of ten times the number of cells 
on one side of the problem space. More precisely for cubical 

cells it takes time steps to traverse a single cell when the 
time step is set by the Courant stability condition 

Ax = cell size dimension (13) 



The condition on N is then that 

1 1 

number cells on a side 
of the problem space 


N 


10 x ( 


/ 3N 


C 3 ) 


NC 


(14) 


The earliest aircraft modeling using FDTD with approximately 30 
cells on a side required approximately 500 time steps. For more 
recent modeling with approximately 100 cells on a side, 2000 or 
more time steps are used. 

For (100 cell) 3 problem spaces, 24 MBytes of memory are 
required to store the fields. Problems on the order of this size 
have been run on a Silicon Graphics 4D 220 with 32 MBytes of 
memory, IBM RISC 6000, an Intel 486 based machine, and VAX 
11/785. Storage is only a problem as in the case of the 486 
where only 16 MBytes of memory was available. This limited the 
problem space size to approximately (80 cells) . 
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For (100 cell) problems with approximately 2000 time steps, 
there is a total of 120 x 10 9 operations to perform. The speeds 
of the previously mentioned machines are 24 MFLOPs (4 processor 
upgraded version), 10 MFLOPS, 1.5 MFLOPS, and 0.2 MFLOPs. T^e 
run times are then 5 3 x 10 3 seconds, 12 x 10 3 seconds, 80 x 10 
seconds and 600 x 10 seconds, respectively. In hours the times 
are 1.4, 3.3, 22.2 and 167 hours. Problems of this size are 
possible on all but the last machine and can in fact be performed 
on a personal computer (486) if one day turnarounds are 
permissible. 

V. VERSION A CODE CAPABILITIES 

The Penn State University FDTD Electromagnetic Scattering 
Code Version A has the following capabilities: 

1) Ability to model lossy dielectric and perfectly conducting 
scatterers. 

2) First and second order outer radiation boundary condition 
(ORBC) operating on the electric fields for dielectric or 
perfectly conducting scatterers. 

3) Near to far zone transformation capability to obtain far zone 
scattered fields. 

4) Gaussian and smooth cosine incident waveforms with arbitrary 
incidence angles. 

5) Near zone field, current or power sampling capability. 

6) Companion code for computing Radar Cross Section (RCS) . 

VI. DEFAULT SCATTERING GEOMETRY 

The code as delivered is set up to calculate the far zone 
backscatter fields for an infinitely thin, 29 cm square, 
perfectly conducting plate. The problem space size is 60 by 60 
by 49 cells in the x, y and z directions, the cells are l cm 
cubes, and the incident waveform is a 0-polarized Gaussian pulse 
with incidence angles of 0=45 and 0=30 degrees. The output data 
files are included as a reference along with a code (RCS3D.F0R) 
for computing the frequency domain RCS using these output data 
files. The ORBC is the second order absorbing boundary condition 
set forth by Mur [2]. 

VII . SUBROUTINE DESCRIPTION 

In the description for each subroutine, an asterisk (*) will 
be placed by the subroutine name if that particular subroutine is 
normally modified when defining a scattering problem. 
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MAIN ROUTINE 

The main routine in the program contains the calls for all 
necessary subroutines to initialize the problem space and 
scattering object (s) and for the incident waveform, far zone 
transformation, field update subroutines, outer radiation 
boundary conditions and field sampling. 

The main routine begins with the include statement and then 
appropriate data files are opened, and subroutines ZERO, BUILD 
and SETUP are called to initialize variables and/or arrays, build 
the object (s) and initialize the incident waveform and 
miscellaneous parameters, respectively. Subroutine SETFZ is 
called to intialize parameters for the near to far zone 
transformation if far zone fields are desired. 

The main loop is entered next, where all of the primary 
field computations and data saving takes place. During each time 
step cycle, the EXSFLD, EYSFLD, and EZSFLD subroutines are called 
to update the x, y, and z components of the scattered electric 
field. The six electric field outer radiation boundary 
conditions (RADE??) are called next to absorb any outgoing 
scattered fields. Time is then advanced 1/2 time step according 
to the Yee algorithm and then the HXSFLD, HYSFLD , AND HZSFLD 
subroutines are called to update the x, y, and z components of 
scattered magnetic field. Time is then advanced another 1/2 step 
and then either near zone fields are sampled and written to disk 
in DATSAV, and/or the near zone to far zone vector potentials are 
updated in SAVFZ. The parameter NZFZ (described later) in the 
common file defines the type of output fields desired. 

After execution of all time steps in the main field update 
loop, subroutine FAROUT is called if far zone fields are desired 
to compute the far zone fields and write them to disk. At this 
point, the execution is complete. 

SUBROUTINE SETFZ 

This subroutine initializes the necessary parameters 
required for far zone field computations. The code as furnished 
computes backscatter far zone fields and can compute bistatic far 
zone fields for one scattering angle (i.e. one 0 and <p angle) . 
Refer to reference [3] for a complete description of the near to 
far zone transformation. Other versions of this subroutine 
provide for multiple bistatic angles. 

SUBROUTINE SAVFZ 

This subroutine updates the near zone to far zone vector 
potentials. 
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SUBROUTINE FAROUT 

This subroutine changes the near zone to far zone vector 
potentials to far zone electric field 0 and <p components and 
writes them to disk. 

SUBROUTINE BUILD * 

This subroutine "builds" the scattering object (s) by 
initializing the I DONE , IDTWO, and IDTHRE arrays. The 
IDONE-IDTHRE arrays are for specifying perfectly conducting and 
lossy dielectric materials. Refer to Figure 1 for a diagram of 
the basic Yee cell. For example, setting an element of the IDONE 
array at some I,J,K location is actually locating dielectric 
material at a cell edge whose center location is I+0.5,J,K. 

Thus, materials with diagonal permittivity tensors can be 
modeled. The default material type for all ID??? arrays is 0, or 
free space. By initializing these arrays to values other than 0, 
the user is defining an object by determining what material types 
are present at each spatial location. Other material types 
available for IDONE-IDTHRE are 1 for perfectly conducting objects 
and 2-9 for lossy non-magnetic dielectrics. It is assumed 
throughout the code that all dielectric materials are non- 
magnetic (i.e. the materials have a permeability of m 0 ) • This 
subroutine also has a section that checks the ID??? arrays to 
determine if legal material types have been defined throughout 
the problem space. The actual material parameters (e and a) are 
defined in subroutine SETUP. The default geometry is a 29 cm 
square perfectly conducting plate. 

The user must be careful that his/her object created in the 
BUILD subroutine is properly formed. 

When it is important to place the object in the center of 
the problem space (to have lowest possible cross-pol scattering 
for symmetric objects) , NX etc. should be odd. This is due to 
the field locations in the Yee cell and also the placement of the 
E field absorbing boundary condition surfaces. 

If the object being modeled has curved surfaces, edges, etc. 
that are at an angle to one or more of the coordinate axes , then 
that shape must be approximately modeled by lines and faces in a 
"stair-stepped" (or stair-cased) fashion. This stair-cased 
approximation introduces errors into computations at higher 
frequencies. Intuitively, the error becomes smaller as more 
cells are used to stair-case a particular object. 

SUBROUTINE DCUBE 

This subroutine builds cubes of dielectric material by 
defining four each of IDONE, IDTWO and IDTHRE components 
corresponding to one spatial cube of dielectric material. It can 
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also be used to define thin (i.e. up to one cell thick) 
dielectric or perfectly conducting plates. Refer to comments 
within DCUBE for a description of the arguments and usage of the 
subroutine. 

SUBROUTINE SETUP * 

This subroutine initializes many of the constants required 
for incident field definition, field update equations, outer 
radiation boundary conditions and material parameters. The 
material parameters e and a are defined for each material type 
using the material arrays EPS and SIGMA respectively. The array 
EPS is used for the total permittivity and SIGMA is used for the 
electric conductivity. These arrays are initialized in SETUP to 
free space material parameters for all material types and then 
the user is required to modify these arrays for his/her 
scattering materials. Thus, for the lossy dielectric material 
type 2, the user must define EPS(2) and SIGMA(2) . The remainder 
of the subroutine computes constants used in field update 
equations and boundary conditions and writes the diagnostics 
file. 

SUBROUTINE EXSFLD 

This subroutine updates all x components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDONE array are used to determine the type of material present 
and the corresponding update equation to be used. These 
scattered field equations are based upon the development given in 
[4]. 

SUBROUTINE EYSFLD 

This subroutine updates all y components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDTWO array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINE EZSFLD 

This subroutine updates all z components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDTHRE array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINES RADEYX , RADEZX, RADEZY, RADEXY, RADEXZ and RADEYZ 

These subroutines apply the outer radiation boundary 
conditions to the scattered electric field on the outer 


boundaries of the problem space. 

SUBROUTINE HXSFLD 

This subroutine updates all x components of scattered 
magnetic field at each time step. The standard non-magnet ic 
update equation is used. 

SUBROUTINE HYSFLD 

This subroutine updates all y components of scattered 
magnetic field at each time step. The standard non-magnetic 
update equation is used. 

SUBROUTINE HZSFLD 

This subroutine updates all z components of scattered 
magnetic field at each time step. The standard non-magnetic 
update equation is used. 

SUBROUTINE DATSAV * 

This subroutine samples near zone scattered field quantities 
and saves them to disk. This subroutine is where the quantities 
to be sampled and their spatial locations are to be specified and 
is only called if near zone f ields only are de sired or if both 
near, an d far zone fields are desired . Total field quantities can 
also be sampled. See comments within the subroutine for 
specifying sampled scattered and/or total field quantities. When 
sampling magnetic fields, remember the <5t/2 time difference 
between E and H when writing the fields to disk. Sections of 
code within this subroutine determine if the sampled quantities 
and the spatial locations have been properly defined. 

FUNCTIONS EXI, EYI and EZI 

These functions are called to compute the x, y and z 
components of incident electric field. The functional form of 
the incident field is contained in a separate function SOURCE. 

FUNCTION SOURCE * 

This function contains the functional form of the incident 
field. The code as furnished uses the Gaussian form of the 
incident field. An incident smooth cosine pulse is also 
available by uncommenting the required lines and commenting out 
the Gaussian pulse. Thus, this function need only be modified if 
the user changes the incident pulse from Gaussian to smooth 
cosine. A slight improvement in computing speed and 
vectorization may be achieved by moving this function inside each 
of the incident field functions EXI, EYI and so on. 
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FUNCTIONS DEXI, DEYI and DEZI 

These functions are called to compute the x, y and z 
components of the time derivative of incident electric field. 

The functional form of the incident field is contained in a 
separate function DSRCE. 

FUNCTION DSRCE * 

This function contains the functional form of the time 
derivative of the incident field. The code as furnished uses the 
time derivative of the Gaussian form of the incident field. A 
smooth cosine pulse time derivative is also available by 
uncommenting the required lines and commenting out the Gaussian 
pulse. Thus, the function need only be modified if the user 
changes from the Gaussian to smooth cosine pulse. Again, a 
slight improvement in computing speed and vectorization may be 
achieved by moving this function inside each of the time 
derivative incident field functions DEXI, DEYI and so on. 

SUBROUTINE ZERO 

This subroutine initializes various arrays and variables to 

zero. 

VIII. INCLUDE PILE DESCRIPTION ( COMMONA. FOR) * 

The include file, COMMONA. FOR, contains all of the arrays 
and variables that are shared among the different subroutines. 
This file will require the most modifications when defining 
scattering problems. A description of the parameters that are 
normally modified follows. 

The parameters NX, NY and NZ specify the size of the problem 
space in cells in the x, y and z directions respectively. For 
problems where it is crucial to center the object within the 
problem space, then NX, NY and NZ should be odd. The parameter 
NTEST defines the number of near zone quantities to be sampled 
and NZFZ defines the field output format. Set NZFZ=0 for near 
zone fields only, NZFZ=1 for far zone fields only and NZFZ=2 for 
both near and far zone fields. Parameter NSTOP defines the 
maximum number of time steps. DELX, DELY, and DELZ (in meters) 
define the cell size in the x, y and z directions respectively. 
The 0 and <p incidence angles (in degrees) are defined by THINC 
and PHINC respectively and the polarization is defined by ETHINC 
and EPHINC. ETHINC=1.0, EPHINC=0 . 0 for 0-polarized incident 
field and ETHINC=0.0, EPHINC=l.O for 0-polarized incident fields. 
Parameters AMP and BETA define the maximum amplitude and the e 
temporal width of the incident pulse respectively. BETA 
automatically adjusts when the cell size is changed and normally 
should not be changed by the user. The far zone scattering 
angles are defined by THETFZ and PHIFZ. The code as furnished 
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performs backscatter computations, but these parameters could be 
modified for a bistatic computation. 

IX. RCS COMPUTATIONS 

A companion code, RCS3D.FOR, has been included to compute 
RCS versus frequency. It uses the file name of the FDTD far zone 
output data (FZ0UT3D.DAT) and writes a data file of far zone 
electric fields versus time (FZTIME.DAT) and RCS versus frequency 
(3DRCS.DAT). The RCS computations are performed up to the 10 
cell/l 0 frequency limit. Refer to comments within this code for 
further details. 

X. RESULTS 

As previously mentioned, the code as furnished models an 
infinitely thin, 29 cm square, perfectly conductinq plate and 
computes backscatter far zone scattered fields at angles of 0=45 
and 0=30 degrees. 

Figures 2-3 shows the co-polarized far zone electric field 
versus time and the co-polarized RCS for the 29 cm square 
perfectly conducting plate. 

Figures 4-5 shows the cross-polarized far zone electric 
field versus time and the cross-polarized RCS for the 29 cm 
square perfectly conducting plate. 

XI. SAMPLE PROBLEM SETUP 

The code as furnished models an infinitely thin, 29 cm 
square, perfectly conducting plate and computes backscatter far 
zone scattered fields at angles of 0=45 and 0=30 degrees. The 
corresponding output data files are also provided, along with a 
code to compute RCS using these data files. In order to change 
the code to a new problem, many different parameters need to be 
modified. A sample problem setup will now be discussed. 

Suppose that the problem to be studied is RCS backscatter 
versus frequency from a 28 cm by 31 cm perfectly conducting plate 
with a 3 cm dielectric coating with a dielectric constant of 4e 0 
using a 0-polarized field. The backscatter angles are 0=30.0 and 
0=60 . 0 degrees and the frequency range is up to 3 Ghz . 

Since the frequency range is up to 3 Ghz, the cell size must 
be chosen appropriately to resolve the field IN ANY MATERIAL at 
the highest frequency of interest. A general rule is that the 
cell size should be 1/10 of the wavelength at the highest 
frequency of interest. For difficult geometries, 1/20 of a 
wavelength may be necessary. The free space wavelength at 3 GHz 
is A. 0 =10 cm and the wavelength in the dielectric coating at 3 GHz 
is 5 cm. The cell size is chosen as 1 cm, which provides a 
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resolution of 5 cells/1 in the dielectric coating and 10 cells/A 0 
in free space. Numerical studies have shown that choosing the 

cell size < 1/4 of the shortest wavelength in any material is 
the practical lower limit. Thus the cell size of 1 cm is barely 
adequate. The cell size in the x, y and z directions is set in 
the common file through variables DELX, DELY and DELZ. Next the 
problem space size must be large enough to accomodate the 
scattering object, plus at least a five cell boundary (10 cells 
is more appropriate) on every side of the object to allow for the 
far zone field integration surface. It is advisable for plate 
scattering to have the plate centered in the x and y directions 
of the problem space in order to reduce the cross-polarized 
backscatter and to position the plate low in the z direction to 
allow strong specular reflections multiple encounters with the 
ORBC . A 10 cell border is chosen, and the problem space size is 
chosen as 49 by 52 by 49 cells in the x, y and z directions 
respectively. As an initial estimate, allow 2048 time steps so 
that energy trapped within the dielectric layer will radiate. 

Thus parameters NX, NY and NZ in COMMONA.FOR would be changed to 
reflect the new problem space size, and parameter NSTOP is 
changed to 2048. If all transients have not been dissipated 
after 2048 time steps, then NSTOP will have to be increased. 
Truncating the time record before all transients have dissipated 
will corrupt frequency domain results. Parameter NZFZ must be 
equal to 1 since we are interested in far zone fields only. To 
build the object, the following lines are inserted into the BUILD 
subroutine: 

C 

C BUILD THE DIELECTRIC SLAB FIRST 

C 

ISTART=11 

JSTART=11 

KSTART=11 

NXWIDE=28 

NYWIDE=31 

NZWIDE=3 

MTYPE=2 

CALL DCUBE ( ISTART , JSTART , KSTART , NXWIDE , NYWIDE , NZWIDE , MTYPE ) 
C 

C BUILD PEC PLATE NEXT 

C 

ISTART=11 

JSTART=11 

KSTART=11 

NXWIDE=28 

NYWIDE=3 1 

NZWIDE=0 

MTYPE=1 

CALL DCUBE (ISTART , JSTART , KSTART , NXWIDE , NYWIDE , NZWIDE , MTYPE ) 
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The PEC plate is built last on the bottom of the dielectric 
slab to avoid any air gaps between the dielectric material and 
the PEC plate. In the common file, the incidence angles THINC 
and PHINC have to be changed to 30.0 and 60.0 respectively, the 
cell sizes (DELX, DELY, DELZ) are set to 0.01, and the 
polarization is set to ETHINC=1.0 and EPHINC=0.0 for 0-polarized 
fields. Since dielectric material 2 is being used for the 
dielectric coating, the constitutive parameters EPS (2) and 
SIGMA ( 2 ) are set to 4e 0 and 0.0 respectively, in subroutine 
SETUP. This completes the code modifications for the sample 
problem. 

XII. NEW PROBLEM CHECKLIST 

This checklist provides a quick reference to determine if 
all parameters have been defined properly for a given scattering 
problem. A reminder when defining quantities within the code: 
use MKS units and specify all angles in degrees. 

COMMONA . FOR : 

1) Is the problem space sized correctly? (NX, NY, NZ) 

2) For near zone fields, is the number of sample points correct? 
(NTEST) 

3) Is parameter NZFZ defined correctly for desired field 
outputs? 

4) Is the number of time steps correct? (NSTOP) 

5) Are the cell dimensions (DELX, DELY, DELZ) defined correctly? 

6) Are the incidence angles (THINC, PHINC) defined correctly? 

7) Is the polarization of the incident wave defined correctly 
(ETHINC, EPHINC)? 

8) For other than backscatter far zone field computations, are 
the scattering angles set correctly? (THETFZ, PHIFZ) 

SUBROUTINE BUILD: 

1) Is the object completely and correctly specified? 

SUBROUTINE SETUP: 

1) Are the constitutive parameters for each material specified 
correctly? (EPS and SIGMA) 
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FUNCTIONS SOURCE and DSRCE: 

1) If the Gaussian pulse is not desired, is it commented out and 

the smooth cosine pulse uncommented? 

SUBROUTINE DATSAV : 

1) For near zone fields, are the sampled field types and spatial 

locations correct for each sampling point? (NTYPE, IOBS, JOBS, 

KOBS) 
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XIV. FIGURE TITLES 

Fig. 1 Standard three dimensional Yee cell showing placement 
of electric and magnetic fields. 

Fig. 2 Co-polarized far zone scattered field versus time for 
29 cm square perfectly conducting plate. 

Fig. 3 Co-polarized RCS versus frequency for 29 cm square 
perfectly conducting plate. 

Fig. 4 Cross-polarized far zone scattered field versus time 
for 29 cm square perfectly conducting plate. 

Fig. 5 Cross-polarized RCS versus frequency for 29 cm square 
perfectly conducting plate. 
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I . INTRODUCTION 

The Penn State Finite Difference Time Domain Electromagnetic 
Scattering Code Version B is a three dimensional numerical 
electromagnetic scattering code based upon the Finite Difference 
Time Domain Technique (FDTD) . The supplied version of the code 
is one version of our current three dimensional FDTD code set. 
This manual provides a description of the code and corresponding 
results for several scattering problems. The manual is organized 
into fourteen sections: introduction, description of the FDTD 

method, operation, resource requirements. Version B code 
capabilities, a brief description of the default scattering 
geometry, a brief description of each subroutine, a description 
of the include file (COMMONB. FOR) , a section briefly discussing 
Radar Cross Section (RCS) computations, a section discussing some 
scattering results, a sample problem setup section, a new problem 
checklist, references and figure titles. 

II . FDTD METHOD 

The Finite Difference Time Domain (FDTD) technique models 
transient electromagnetic scattering and interactions with 
objects of arbitrary shape and/or material composition. The 
technique was first proposed by Yee [1] for isotropic, non- 
dispersive materials in 1966; and has matured within the past 
twenty years into a robust and efficient computational method. 

The present FDTD technique is capabable of transient 
electromagnetic interactions with objects of arbitrary and 
complicated geometrical shape and material composition over a 
large band of frequencies. This technique has recently been 
extended to include dispersive dielectric materials, chiral 
materials and plasmas. 

In the FDTD method, Maxwell's curl equations are discretized 
in time and space and all derivatives (temporal and spatial) are 
approximated by central differences. The electric and magnetic 
fields are interleaved in space and time and are updated in a 
second-order accurate leapfrog scheme. The computational space 
is divided into cells with the electric fields located on the 
edges and the magnetic fields on the faces (see Figure 1). FDTD 
objects are defined by specifying dielectric and/or magnetic 
material parameters at electric and/or magnetic field locations. 

Two basic implementations of the FDTD method are widely used 
for electromagnetic analysis: total field formalism and scattered 
field formalism. In the total field formalism, the electric and 
magnetic field are updated based upon the material type present 
at each spatial location. In the scattered field formalism, the 
incident waveform is defined analytically and the scattered field 
is coupled to the incident field through the different material 
types. For the incident field, any waveform, angle of incidence 
and polarization is possible. The separation of the incident and 
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scattered fields conveniently allows an absorbing boundary to be 
employed at the extremities of the discretized problem space to 
absorb the scattered fields. 

This code is a scattered field code, and the total E and H 
fields may be found by combining the incident and scattered 
fields. Any type of field quantity (incident, scattered, or 
total) , Poynting vector or current are available anywhere within 
the computational space. These fields, incident, scattered and 
total, may be found within, on or about the interaction object 
placed in the problem space. By using a near to far field 
transformation, far fields can be determined from the near fields 
within the problem space thereby affording radiation patterns and 
RCS values. The accuracy of these calculations is typically 
within a dB of analytic solutions for dielectric and magnetic 
sphere scattering. Further improvements are expected as better 
absorbing boundary conditions are developed and incorporated. 


III. OPERATION 


Typically, a truncated Gaussian incident waveform is used to 
excite the system being modeled, however certain code versions 
also provide a smooth cosine waveform for convenience in modeling 
dispersive materials. The interaction object is defined in the 
discretized problem space with arrays at each cell location 
created by the discretization. All three dielectric material 
types for E field components within a cell can be individually 
specified by the arrays IDONE(I, J,K) , IDTWO(I, J,K) , 

IDTHRE (I , J, K) . This models arbitrary dielectric materials with m 
= n 0 . By an obvious extension to six arrays, magnetic materials 
with fi f n Q can be modeled. 

Scattering occurs when the incident wave, marched forward in 
time in small steps set by the Courant stability condition, 
reaches the interaction object. Here a scattered wave must 
appear along with the incident wave so that the Maxwell equations 
are satisfied. If the material is a perfectly conductive metal 
then only the well known boundary condition 


scat inc 

^tan = - ^*tan 


( 1 ) 


must be satisfied. For a nondispersive dielectric the 
requirement is that the total field must satisfy the Maxwell 
equations in the material: 


Vx E tot 


Vx (E inc + E scat 


1 dH tot 
dt 


i a(H inc +H scat ) 


at 


(2) 
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VxH tot = Vx (H inc +H scat ) 


aE tot 

= e 

at 


+ 0 E 


tot 


(3) 


= e a (E inc +E sc ^2 +g(E inc +E scat ) 

at 


(4) 


Additionally the incident wave, defined as moving unimpeded 
through a vacuum in the problem space, satisfies everywhere in 
the problem the Maxwell equations for free space 


Vx E inc 


1 aH inc 
3 t 


(5) 


VxH inc -t 


aE 


inc 


at 


( 6 ) 


Subtracting the second set of equations from the first yields the 
Maxwell equations governing the scattered fields in the material. 


Vx E scat 


1 aH scat 
m 0 at 


(7) 


VxH scat = (e -e ) 


aE’ 


at 


+ aE ,nc +e 


^t? SC8t 

Or, _,scat 

+ CXE 

at 


( 8 ) 


Outside the material this simplifies to: 


„ scat 1 aff 

Vx E 


scat 


M 0 3t 


( 9 ) 


VxH scat =e 


aE ! 


,scat 


at 


( 10 ) 


Magnetic materials, dispersive effects, non-linearities, 
etc., are further generalizations of the above approach. Based 
on the value of the material type, the subroutines for 
calculating scattered E and H field components branch to the 
appropriate expression for that scattered field component and 
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that component is advanced in time according to the selected 
algorithm. As many materials can be modeled as desired, the 
number equals the dimension selected for the flags. If materials 
with behavior different from those described above must be 
modeled, then after the appropriate algorithm is found, the 
code's branching structure allows easy incorporation of the new 
behavior. 

IV. RESOURCE REQUIREMENTS 

The number of cells the problem space is divided into times 
the six components per cell set the problem space storage 
requirements 

Storage = NC x 6 components/ cell x 4 bytes /component (11) 


and the computational cost 

Operations = NC x 6 comp/cell x 10 ops/component x N (12) 


where N is the number of time steps desired. 

N typically is on the order of ten times the number of cells 
on one side of the problem space. More precisely for cubical 

cells it takes time steps to traverse a single cell when the 
time step is set by the Courant stability condition 



Ax = cell size dimension 


(13) 


The condition on N is then that 

1 

3 number cells on a side (14) 

of the problem space 

The earliest aircraft modeling using FDTD with approximately 30 
cells on a side required approximately 500 time steps. For more 
recent modeling with approximately 100 cells on a side, 2000 or 
more time steps are used. 

For (100 cell) 3 problem spaces, 24 MBytes of memory are 
required to store the fields. Problems on the order of this size 
have been run on a Silicon Graphics 4D 220 with 32 MBytes of 
memory, IBM RISC 6000, an Intel 486 based machine, and VAX 
11/785. Storage is only a problem as in the case of the 486 
where only 16 MBytes of memory was available. This limited the 
problem space size to approximately (80 cells) . 


ix ( /7i 


NC 3 ) 
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For (100 cell) 3 problems with approximately 2000 time steps, 
there is a total of 120 x 10 operations to perform. The speeds 
of the previously mentioned machines are 24 MFLOPs (4 processor 
upgraded version), 10 MFLOPS, 1.5 MFLOPS, ^nd 0.2 MFLOPs. T^e 
run times are then 5 x 10 seconds, 12 x 10 seconds, 80 x 10 
seconds and 600 x 10 3 seconds, respectively. In hours the times 
are 1.4, 3.3, 22.2 and 167 hours. Problems of this size are 
possible on all but the last machine and can in fact be performed 
on a personal computer (486) if one day turnarounds are 
permissible. 

V. VERSION B CODE CAPABILITIES 

The Penn State University FDTD Electromagnetic Scattering 
Code Version B has the following capabilities: 

1) Ability to model lossy dielectric and perfectly conducting 
scatterers . 

2) Ability to model dispersive dielectric scatterers. This 
dispersive FDTD method is now designated (FD) TD for Frequency- 
Dependent Finite Difference Time Domain. 

3) First and second order outer radiation boundary condition 
( ORBC ) operating on the electric fields. 

4) Near to far zone transformation capability to obtain far zone 
scattered fields. 

5) Gaussian and smooth cosine incident waveforms with arbitrary 
incidence angles. 

6) Near zone field, current or power sampling capability. 

7) Companion code for computing Radar Cross Section (RCS) . 

VI. DEFAULT SCATTERING GEOMETRY 

The code as delivered is set up to calculate the far zone 
backscatter fields for a 20 cm radius, dispersive, 0.25 dB 
dielectric foam sphere. The 0.25 dB loaded foam is defined by a 
frequency dependent permittivity with an effective DC 
conductivity given by 


e(o) € s -€ « ° 

— — = + + 

€ 0 1+jWTo j W6 0 


(15) 


where e 0 is the infinite frequency permittivity, e s is the zero 
frequency permittivity, r Q is the relaxation time and to is the 
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radian frequency. The 0.25 dB loaded foam parameters are 
e =1.06, 6=1.16, T =0.6497 ns and a=2.954e-4. The problem space 

*0 1 S f U 

size is 65 by 65 by 65 cells in the x, y and z directions, the 
cells are 1 cm cubes, and the incident waveform is a 0-polarized 
smooth cosine pulse with incidence angles of 0=22.5 and 0=22.5 
degrees. The output data files are included as a reference along 
with a code (RCS3D. FOR) for computing the frequency domain RCS 
using these output data files. The ORBC is the second order 
absorbing boundary condition set forth by Mur [2]. 

VII. SUBROUTINE DESCRIPTION 

In the description for each subroutine, an asterisk (*) will 
be placed by the subroutine name if that particular subroutine is 
normally modified when defining a scattering problem. 

MAIN ROUTINE 

The main routine in the program contains the calls for all 
necessary subroutines to initialize the problem space and 
scattering object (s) and for the incident waveform, far zone 
transformation, field update subroutines, outer radiation 
boundary conditions and field sampling. 

The main routine begins with the include statement and then 
appropriate data files are opened, and subroutines ZERO, BUILD 
and SETUP are called to initialize variables and/or arrays, build 
the object (s) and initialize the incident waveform and 
miscellaneous parameters, respectively. Subroutine SETFZ is 
called to intialize parameters for the near to far zone 
transformation if far zone fields are desired. 

The main loop is entered next, where all of the primary 
field computations and data saving takes place. During each time 
step cycle, the EXSFLD , EYSFLD, and EZSFLD subroutines are called 
to update the x, y, and z components of the scattered electric 
field. The six electric field outer radiation boundary 
conditions (RADE??) are called next to absorb any outgoing 
scattered fields. Time is then advanced 1/2 time step according 
to the Yee algorithm and then the HXSFLD, HYSFLD, AND HZSFLD 
subroutines are called to update the x, y, and z components of 
scattered magnetic field. Time is then advanced another 1/2 step 
and then either near zone fields are sampled and written to disk 
in DATSAV, and/or the near zone to far zone vector potentials are 
updated in SAVFZ. The parameter NZFZ (described later) in the 
common file defines the type of output fields desired. 

After execution of all time steps in the main field update 
loop, subroutine FAROUT is called if far zone fields are desired 
to compute the far zone fields and write them to disk. At this 
point, the execution is complete. 
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SUBROUTINE SETFZ 

This subroutine initializes the necessary parameters 
required for far zone field computations. The code as furnished 
computes backscatter far zone fields and can compute bistatic far 
zone fields for one scattering angle (i.e. one 0 and <p angle) . 
Refer to reference [3] for a complete description of the near to 
far zone transformation. Other versions of this subroutine 
provide for multiple bistatic angles. 

SUBROUTINE SAVFZ 

This subroutine updates the near zone to far zone vector 
potentials. 

SUBROUTINE FAROUT 

This subroutine changes the near zone to far zone vector 
potentials to far zone electric field 0 and (p components and 
writes them to disk. 

SUBROUTINE BUILD * 

This subroutine "builds" the scattering object (s) by 
initializing the I DONE, IDTWO, and IDTHRE arrays to specify 
perfectly conducting, lossy dielectrics and dispersive dielectric 
materials. Refer to Figure 1 for a diagram of the basic Yee 
cell. Setting an element of the IDONE array at some I,J,K 
location is actually locating dielectric material at a cell edge 
whose center location is I+0.5,J,K. Thus, materials with 
diagonal permittivity tensors can be modeled. The default 
material type for all ID??? arrays is 0, or free space. By 
initializing these arrays to values other than 0, the user is 
defining an object by determining what material types are present 
at each spatial location. Other material types available for 
IDONE- IDTHRE are 1 for perfectly conducting objects, 2-9 for 
lossy non-magnetic dielectrics, 20-29 for dispersive dielectrics. 
This code assumes that all dielectric and dispersive dielectric 
materials have a permeability of n Q . This subroutine also has a 
section that checks the ID??? arrays to determine if legal 
material types have been defined throughout the problem space. 

The actual non-dispersive material parameters (e, jli, and a) are 
defined in subroutine SETUP. The dispersive material parameters 

(e s , e ffl , r Q , a) are also defined in a separate section in SETUP. 

The default geometry is a 20 cm radius, dispersive, 0.25 dB 
dielectric foam sphere. 

The user must be careful that his/her object created in the 
BUILD subroutine is properly formed. 
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When it is important to place the object in the center of 
the problem space (to have lowest possible cross-pol scattering 
for symmetric objects) , NX etc. should be odd. This is due to 
the field locations in the Yee cell and also the placement of the 
E field absorbing boundary condition surfaces. 

If the object being modeled has curved surfaces, edges, etc. 
that are at an angle to one or more of the coordinate axes, then 
that shape must be approximately modeled by lines and faces in a 
’’stair— stepped" (or stair— cased) fashion. This stair— cased 
approximation introduces errors into computations at higher 
frequencies. Intuitively, the error becomes smaller as more 
cells are used to stair-case a particular object. Thus, the 
default 0.25 dB dielectric foam sphere scattering geometry is a 
stair-cased sphere. 

SUBROUTINE DCUBE 

This subroutine builds cubes of dielectric material by 
defining four each of I DONE, IDTWO and IDTHRE components 
corresponding to one spatial cube of dielectric material. It can 
also be used to define thin (i.e. up to one cell thick) 
dielectric or perfectly conducting plates. Refer to comments 
within DCUBE for a description of the arguments and usage of the 
subroutine. 

SUBROUTINE SETUP * 

This subroutine initializes many of the constants required 
for incident field definition, field update equations, outer 
radiation boundary conditions and material parameters. The 
material parameters e and a are defined for each material type 
(non-dispersive) using the material arrays EPS and SIGMA 
respectively. The array EPS is used for the total permittivity, 
and SIGMA is used for the electric conductivity. These arrays 
are initialized in SETUP to free space material parameters for 
all material types and then the user is required to modify these 
arrays for his/her scattering materials. Thus, for the lossy 
dielectric material type 2, the user must define EPS (2) and 
SIGMA ( 2 ) . 

For dispersive dielectric materials, different material 
parameter arrays are used. The functional form of the frequency 
dependent permittivity that was implemented in the code is the 
Debye relaxation [4] with an effective DC conductivity given by 

£(<*>) = € / - je // = e B e + e 0 X e (<*>) + ~~ 

jto 


( 16 ) 
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where the frequency dependent electric susceptibility function is 
defined as 


x*(«) 


e g - e, 

1 + jux 0 


(17) 


where e m is the infinite frequency permittivity, e s is the zero 

frequency permittivity, r Q is the relaxation time, o is the 

effective electric conductivity, and o> is the radian frequency. 
The corresponding time domain susceptibility function is given by 


X„(t) = 


/£ s “ 6 ^ 




u(t) 


(18) 


The FDTD implementation of frequency dependent permittivity 
involves a convolution with the electric field and interested 
readers are referred to references [5-6] for further details. 

For dispersive dielectric materials, the corresponding 
material parameter arrays are EPSINF (£„), EPSSTA (e s ) , RELAXT 

(T 0 ), and RELSIG (a). These dispersive material parameters are 

defined under the DISPERSIVE SETUP portion of the subroutine. 

The remainder of the subroutine computes constants used in field 
update equations and boundary conditions and writes the 
diagnostics file. 

SUBROUTINE EXSFLD 

This subroutine updates all x components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDONE array are used to determine the type of material present 
and the corresponding update equation to be used. These 
scattered field equations are based upon the development given in 
[7] . 

SUBROUTINE EYSFLD 

This subroutine updates all y components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDTWO array are used to determine the type of material present 
and the corresponding update equation to be used. 
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SUBROUTINE EZSFLD 

This subroutine updates all z components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDTHRE array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINES RADEYX , RADEZX, RADEZY, RADEXY, RADEXZ and RADEYZ 

These subroutines apply the outer radiation boundary 
conditions to the scattered electric field on the outer 
boundaries of the problem space. 

SUBROUTINE HXSFLD 

This subroutine updates all x components of scattered 
magnetic field at each time step. The standard non-magnet ic 
update equation is used. 

SUBROUTINE HYSFLD 

This subroutine updates all y components of scattered 
magnetic field at each time step. The standard non-magnet ic 
update equation is used. 

SUBROUTINE HZSFLD 

This subroutine updates all z components of scattered 
magnetic field at each time step. The standard non-magnetic 
update equation is used. 

SUBROUTINE DATSAV * 

This subroutine samples near zone scattered field quantities 
and saves them to disk. This subroutine is where the quantities 
to be sampled and their spatial locations are to be specified and 
is only called if near zone fields only are desired or if both 
near and far zone fields are desired. Total field quantities can 
also be sampled. See comments within the subroutine for 
specifying sampled scattered and/or total field quantities. When 
sampling magnetic fields, remember the S t/2 time difference 
between E and H when writing the fields to disk. Sections of 
code within this subroutine determine if the sampled quantities 
and the spatial locations have been properly defined. 

FUNCTIONS EXI, EYI and EZI 


These functions are called to compute the x, y and z 
components of incident electric field. The functional form of 
the incident field is contained in a separate function SOURCE. 
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FUNCTION SOURCE * 

This function contains the functional form of the incident 
field. The code as furnished uses the smooth cosine form of the 
incident field. An incident Gaussian pulse is also available by 
uncommenting the required lines and commenting out the smooth 
cosine pulse. Thus, this function need only be modified if the 
user changes the incident pulse from smooth cosine to Gaussian. 
Currently, only the smooth cosine pulse can be used for 
scattering from dispersive targets. A slight improvement in 
computing speed and vectorization may be achieved by moving this 
function inside each of the incident field functions EXI, EYI and 
so on. 

FUNCTIONS DEXI, DEYI and DEZI 

These functions are called to compute the x, y and z 
components of the time derivative of incident electric field. 

The functional form of the incident field is contained in a 
separate function DSRCE. 

FUNCTIONS DEXIXE, DEYIXE and DEZIXE 

These functions compute the x, y and z components of the 
convolution of the time derivative of the incident field with the 

electric Debye susceptibility function % . 

FUNCTION DSRCE * 

This function contains the functional form of the time 
derivative of the incident field. The code as furnished uses the 
time derivative of the smooth cosine form of the incident field. 

A Gaussian pulse time derivative is also available by 
uncommenting the required lines and commenting out the smooth 
cosine pulse. Thus, the function need only be modified if the 
user changes from the smooth cosine to Gaussian pulse. Again, a 
slight improvement in computing speed and vectorization may be 
achieved by moving this function inside each of the time 
derivative incident field functions DEXI, DEYI and so on. 

FUNCTION DCONV 

This function evaluates the convolution of the time 
derivative of the incident field with the Debye susceptibility 

function x • 

SUBROUTINE ZERO 

This subroutine initializes various arrays and variables to 


zero . 
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VIII. INCLUDE FILE DESCRIPTION ( COMMONB . FOR) * 

The include file, COMMONB. FOR, contains all of the arrays 
and variables that are shared among the different subroutines. 
This file will require the most modifications when defining 
scattering problems. A description of the parameters that are 
normally modified follows. 

The parameters NX, NY and NZ specify the size of the problem 
space in cells in the x, y and z directions respectively. For 
problems where it is crucial to center the object within the 
problem space, then NX, NY and NZ should be odd. The parameter 
NTEST defines the number of near zone quantities to be sampled 
and NZFZ defines the field output format. Set NZFZ=0 for near 
zone fields only, NZFZ=1 for far zone fields only and NZFZ=2 for 
both near and far zone fields. Parameter NUMMAT defines the 
total number of material types that are available for use. 

NEDISP defines the number of dispersive dielectric materials that 
are being used. Parameter NSTOP defines the maximum number of 
time steps. DELX, DELY, and DELZ (in meters) define the cell 
size in the x, y and z directions respectively. The 0 and <p 
incidence angles (in degrees) are defined by THINC and PHINC 
respectively and the polarization is defined by ETHINC and 
EPHINC. ETHINC=1.0, EPHINC=0.0 for 0-polarized incident field 
and ETHINC=0.0, EPHINC=1 . 0 for 0-polarized incident fields. 2 
Parameters AMP and BETA define the maximum amplitude and the e 
temporal width of the incident pulse respectively. BETA 
automatically adjusts when the cell size is changed and normally 
should not be changed by the user. The far zone scattering 
angles are defined by THETFZ and PHIFZ. The code as furnished 
performs backscatter computations, but these parameters could be 
modified for a bistatic computation. 

IX. RCS COMPUTATIONS 

A companion code, RCS3D. FOR, has been included to compute 
RCS versus frequency. It uses the file name of the (FD) TD far 
zone output data (FZ0UT3D.DAT) and writes data files of far zone 
electric fields versus time (FZTIME.DAT) and RCS versus frequency 
(3DRCS.DAT). The RCS computations are performed up to the 10 
cell/A 0 frequency limit. Refer to comments within this code for 
further details. 

X. RESULTS 

As previously mentioned, the code as furnished models a 20 
cm radius, dispersive, 0.25 dB dieletric foam sphere and computes 
backscatter far zone scattered fields at angles of 0=22.5 and 
0=22.5 degrees. Results are included for the 0.25 dB dielectric 
foam sphere and for a 60 dB dielectric foam sphere. The material 
parameters for 0.25 dB and 60 dB loaded foam are: 
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0.25 DB FOAM 


60 DB FOAM 


1.16 e s 

1.01 e. 

0.6497 ns r Q 

2 . 954E-04 S/m a 


41.0 

1.6 

0.3450 ns 
3 . 902E-01 S/m 


For the 0.25 dB foam there are 10 cells per wavelength at 
approximately 3.0 GHz. For the 60 dB foam there are 10 cells per 
wavelength at approximately 2.35 GHz, as the 60 dB foam has a 
higher dielectric constant. 


Figures 2-7 show the real and imaginary parts of the 0.25 dB 
and 60 dB foam permeability and the dielectric susceptibilities 
versus time for both materials. 


Figures 8-11 show the co-polarized far zone electric field 
versus time and the co-polarized RCS for the 0.25 dB and 60 dB 
dielectric foam spheres respectively. 


XI. SAMPLE PROBLEM SETUP 

The code as furnished models a 20 cm radius, dispersive, 

0.25 dB dielectric foam sphere and computes backscatter far zone 
scattered fields at angles of 0=22.5 and <£=22.5 degrees. The 
corresponding output data files are also provided, along with a 
code to compute Radar Cross Section using these data files. In 
order to change the code to a new problem, many different 
parameters need to be modified. A sample problem setup will now 
be discussed. 

Suppose that the problem to be studied is RCS backscatter 
versus frequency from a 28 cm by 31 cm perfectly conducting plate 
with a 3 cm dielectric coating with a dielectric constant of 4c 0 
using a 0-polarized field. The backscatter angles are 0=30.0 and 
0=60.0 degrees and the frequency range is up to 3 Ghz. 

Since the frequency range is up to 3 Ghz, the cell size must 
be chosen appropriately to resolve the field IN ANY MATERIAL at 
the highest frequency of interest. A general rule is that the 
cell size should be 1/10 of the wavelength at the highest 
frequency of interest. For difficult geometries, 1/20 of a 
wavelength may be necessary. The free space wavelength at 3 GHz 
is A, 0 =10 cm and the wavelength in the dielectric coating at 3 GHz 
is 5 cm. The cell size is chosen as 1 cm, which provides a 
resolution of 5 cells/A, in the dielectric coating and 10 cells/X 0 
in free space. Numerical studies have shown that choosing the 

cell size < 1/4 of the shortest wavelength in any material is 
the practical lower limit. Thus the cell size of 1 cm is barely 
adequate. The cell size in the x, y and z directions is set in 
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the common file through variables DELX, DELY and DELZ. Next the 
problem space size must be large enough to accomodate the 
scattering object, plus at least a five cell boundary (10 cells 
is more appropriate) on every side of the object to allow for the 
far zone field integration surface. It is advisable for plate 
scattering to have the plate centered in the x and y directions 
of the problem space in order to reduce the cross-polarized 
backscatter and to position the plate low in the z direction to 
allow strong specular reflections multiple encounters with the 
ORBC. A 10 cell border is chosen, and the problem space size is 
chosen as 49 by 52 by 49 cells in the x, y and z directions 
respectively. As an initial estimate, allow 2048 time steps so 
that energy trapped within the dielectric layer will radiate. 

Thus parameters NX, NY and NZ in COMMONB . FOR would be changed to 
reflect the new problem space size, and parameter NSTOP is 
changed to 2048. If all transients have not been dissipated 
after 2048 time steps, then NSTOP will have to be increased. 
Truncating the time record before all transients have dissipated 
will corrupt frequency domain results. Parameter NZFZ must be 
equal to 1 since we are interested in far zone fields only. To 
build the object, the following lines are inserted into the BUILD 
subroutine: 

C 

C BUILD THE DIELECTRIC SLAB FIRST 

C 

ISTART=1 1 

JSTART=11 

KSTART= 1 1 

NXWIDE=28 

NYWIDE=31 

NZWIDE=3 

MTYPE=2 

CALL DCUBE ( I START , JSTART , KSTART , NXWIDE , NYWIDE , NZWIDE , MTYPE) 
C 

C BUILD PEC PLATE NEXT 

C 

ISTART=11 

JSTART=11 

KSTART=11 

NXWIDE=2 8 

NYWIDE=3 1 

NZWIDE=0 

MTYPE=1 

CALL DCUBE ( ISTART , JSTART , KSTART, NXWIDE , NYWIDE , NZWIDE , MTYPE) 


The PEC plate is built last on the bottom of the dielectric 
slab to avoid any air gaps between the dielectric material and 
the PEC plate. In the common file, the incidence angles THINC 
and PHINC have to be changed to 30.0 and 60.0 respectively, the 
cell sizes (DELX, DELY, DELZ) are set to 0.01, and the 
polarization is set to ETHINC=1.0 and EPHINC=0.0 for 0-polarized 
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fields. Since dielectric material 2 is being used for the 
dielectric coating, the constitutive parameters EPS (2) and 
SIGMA ( 2 ) are set to 4e 0 and 0.0 respectively, in subroutine 
SETUP. This completes the code modifications for the sample 
problem. 

XII. NEW PROBLEM CHECKLIST 

This checklist provides a quick reference to determine if 
all parameters have been defined properly for a given scattering 
problem. A reminder when defining quantities within the code: 
use MKS units and specify all angles in degrees. 

COMMONB . FOR : 

1) Is the problem space sized correctly? (NX, NY, NZ) 

2) For near zone fields, is the number of sample points correct? 
(NTEST) 

3) Is parameter NZFZ defined correctly for desired field 
outputs? 

4) Is the number of dispersive dielectric (NEDISP) materials 
defined correctly? 

5) Is the number of time steps correct? (NSTOP) 

6) Are the cell dimensions ( DELX , DELY, DELZ) defined correctly? 

7) Are the incidence angles (THINC, PHINC) defined correctly? 

8) Is the polarization of the incident wave defined correctly 
(ETHINC, EPHINC)? 

9) For other than backscatter far zone field computations, are 
the scattering angles set correctly? (THETFZ, PHIFZ) 

SUBROUTINE BUILD: 

1) Is the object completely and correctly specified? 

SUBROUTINE SETUP: 

1) Are the constitutive parameters for each material specified 
correctly? (EPS and SIGMA) 

2) Are the constitutive parameters for each dispersive material 
defined correctly? (EPSSTA, EPSINF, RELAXT , RELSIG) 
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FUNCTIONS SOURCE and DSRCE: 

1) If the smooth cosine pulse is not desired, is it commented 

out and the Gaussian pulse uncommented? 

SUBROUTINE DATSAV : 

1) For near zone fields, are the sampled field types and spatial 

locations correct for each sampling point? (NTYPE, IOBS, JOBS, 

KOBS) 
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XIV. FIGURE TITLES 

Fig. 1 Standard three dimensional Yee cell showing placement 
of electric and magnetic fields. 

Fig. 2 Real part of relative permittivity versus frequency for 
0.25 dB dielectric foam. 
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Imaginary part of relative permittivity versus 
frequency for 0.25 dB dielectric foam. 

Relative dielectric susceptibility versus time for 0.25 
dB dielectric foam. 

Real part of relative permittivity versus frequency for 
60 dB dielectric foam. 

Imaginary part of relative permittivity versus 
frequency for 60 dB dielectric foam. 

Relative electric susceptibility versus time for 60 dB 
dielectric foam. 

Co-polarized far zone scattered field versus time for 
0.25 dB dielectric foam sphere with 20 cm radius. 

Co-polarized Radar Cross Section versus frequency for 
0.25 dB dielectric foam sphere with 20 cm radius. 

Co-polarized far zone scattered field versus time for 
60 dB dielectric foam sphere with 20 cm radius. 

Co-polarized Radar Cross Section versus frequency for 
60 dB dielectric foam sphere with 20 cm radius. 
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I . INTRODUCTION 

The Penn State Finite Difference Time Domain Electromagnetic 
Scattering Code Version C is a three dimensional numerical 
electromagnetic scattering code based upon the Finite Difference 
Time Domain Technique (FDTD) . The supplied version of the code 
is one version of our current three dimensional FDTD code set. 
This manual provides a description of the code and corresponding 
results for several scattering problems. The manual is organized 
into fourteen sections: introduction, description of the FDTD 

method, operation, resource requirements. Version C code 
capabilities, a brief description of the default scattering 
geometry, a brief description of each subroutine, a description 
of the include file (COMMONC. FOR) , a section briefly discussing 
Radar Cross Section (RCS) computations, a section discussing some 
scattering results, a sample problem setup section, a new problem 
checklist, references and figure titles. 

II. FDTD METHOD 

The Finite Difference Time Domain (FDTD) technique models 
transient electromagnetic scattering and interactions with 
objects of arbitrary shape and/or material composition. The 
technique was first proposed by Yee [1] for isotropic, non- 
dispersive materials in 1966; and has matured within the past 
twenty years into a robust and efficient computational method. 

The present FDTD technique is capabable of transient 
electromagnetic interactions with objects of arbitrary and 
complicated geometrical shape and material composition over a 
large band of frequencies. This technique has recently been 
extended to include dispersive dielectric materials, chiral 
materials and plasmas. 

In the FDTD method, Maxwell's curl equations are discretized 
in time and space and all derivatives (temporal and spatial) are 
approximated by central differences. The electric and magnetic 
fields are interleaved in space and time and are updated in a 
second-order accurate leapfrog scheme. The computational space 
is divided into cells with the electric fields located on the 
edges and the magnetic fields on the faces (see Figure 1) . FDTD 
objects are defined by specifying dielectric and/or magnetic 
material parameters at electric and/or magnetic field locations. 

Two basic implementations of the FDTD method are widely used 
for electromagnetic analysis: total field formalism and scattered 
field formalism. In the total field formalism, the electric and 
magnetic field are updated based upon the material type present 
at each spatial location. In the scattered field formalism, the 
incident waveform is defined analytically and the scattered field 
is coupled to the incident field through the different material 
types. For the incident field, any waveform, angle of incidence 
and polarization is possible. The separation of the incident and 
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scattered fields conveniently allows an absorbing boundary to be 
employed at the extremities of the discretized problem space to 
absorb the scattered fields. 

This code is a scattered field code, and the total E and H 
fields may be found by combining the incident and scattered 
fields. Any type of field quantity (incident, scattered, or 
total) , Poynting vector or current are available anywhere within 
the computational space. These fields, incident, scattered and 
total, may be found within, on or about the interaction object 
placed in the problem space. By using a near to far field 
transformation, far fields can be determined from the near fields 
within the problem space thereby affording radiation patterns and 
RCS values. The accuracy of these calculations is typically 
within a dB of analytic solutions for dielectric and magnetic 
sphere scattering. Further improvements are expected as better 
absorbing boundary conditions are developed and incorporated. 

III. OPERATION 


Typically, a truncated Gaussian incident waveform is used to 
excite the system being modeled, however certain code versions 
also provide a smooth cosine waveform for convenience in modeling 
dispersive materials. The interaction object is defined in the 
discretized problem space with arrays at each cell location 
created by the discretization. All three dielectric material 
types for E field components within a cell can be individually 
specified by the arrays IDONE(I, J,K) , IDTWO(I, J,K) , 

IDTHRE ( I , J , K ) . This models arbitrary dielectric materials with /i 
= Mo* B Y an obvious extension to six arrays, magnetic materials 
with /i f /i 0 can be modeled. 


Scattering occurs when the incident wave, marched forward in 
time in small steps set by the Courant stability condition, 
reaches the interaction object. Here a scattered wave must 
appear along with the incident wave so that the Maxwell equations 
are satisfied. If the material is a perfectly conductive metal 
then only the well known boundary condition 


scat inc 

^tan “ ~^tan 


(i) 


must be satisfied. For a nondispersive dielectric the 
requirement is that the total field must satisfy the Maxwell 
equations in the material: 


Vx E tot 


Vx ( E inc + E scat 


1 8H tot 
M 0 dt 


1 d (H ,nc +H scat ) 

il 


at 


( 2 ) 
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VxH tot = Vx (H ,nc +H scat ) = € 


inc , TT scatx 


3E 


tot 


dt 


+ o E 


tot 


(3) 


= e 


<9 ( E inc + E scat ) 

at 


+ a (E inc + E scat ) 


(4) 


Additionally the incident wave, defined as moving unimpeded 
through a vacuum in the problem space, satisfies everywhere in 
the problem the Maxwell equations for free space 


Vx E inc 


i aH inc 
M 0 dt 


(5) 


VxH inc 


= e 


aE inc 

0 at 


( 6 ) 


Subtracting the second set of equations from the first yields the 
Maxwell equations governing the scattered fields in the material: 


Vx E scat 


1 aH scat 
M 0 dt 


(7) 


aE' 


aE 


scat 


VxH = (e -e ) — — +aE""‘ +c 

0 at at 


+ crE 


( 8 ) 


Outside the material this simplifies to: 


_ craf 1 aH 

Vx E scat = - — 


scat 


M 0 dt 


( 9 ) 


Vx H scat =e 


aE 


scat 


at 


( 10 ) 


Magnetic materials, dispersive effects, non-linearities, 
etc. , are further generalizations of the above approach. Based 
on the value of the material type, the subroutines for 
calculating scattered E and H field components branch to the 
appropriate expression for that scattered field component and 


7 


that component is advanced in time according to the selected 
algorithm. As many materials can be modeled as desired, the 
number equals the dimension selected for the flags. If materials 
with behavior different from those described above must be 
modeled, then after the appropriate algorithm is found, the 
code's branching structure allows easy incorporation of the new 
behavior. 

IV. RESOURCE REQUIREMENTS 

The number of cells the problem space is divided into times 
the six components per cell set the problem space storage 
requirements 

Storage = NC x 6 components/cell x 4 bytes /component (H) 

and the computational cost 

Operations = NC x 6 comp/cell x io ops/component x N (12) 

where N is the number of time steps desired. 

N typically is on the order of ten times the number of cells 
on one side of the problem space. More precisely for cubical 

cells it takes >/3 

time steps to traverse a single cell when the 
time step is set by the Courant stability condition 

Ax 

At = Ax = cell size dimension (13) 

y 3 c 


The condition on N is then that 

1 1 

number cells on a side 
of the problem space 


N - 10 x ( 


/ 7 ; 


NC 3 ) 


NC 


(14) 


The earliest aircraft modeling using FDTD with approximately 30 
cells on a side required approximately 500 time steps. For more 
recent modeling with approximately 100 cells on a side, 2000 or 
more time steps are used. 

For (100 cell) 3 problem spaces, 24 MBytes of memory are 
required to store the fields. Problems on the order of this size 
have been run on a Silicon Graphics 4D 220 with 32 MBytes of 
memory, IBM RISC 6000, an Intel 486 based machine, and VAX 
11/785. Storage is only a problem as in the case of the 486 
where only 16 MBytes of memory was available. 3 This limited the 
problem space size to approximately (80 cells) . 



8 


For (100 cell) problem^ with approximately 2000 time steps, 
there is a total of 120 x 10 v operations to perform. The speeds 
of the previously mentioned machines are 24 MFLOPs (4 processor 
upgraded version), 10 MFLOPS, 1.5 MFLOPS, and 0.2 MFLOPs. The 
run times are then 5 3 x 10 3 seconds, 12 x 10 3 seconds, 80 x 10 3 
seconds and 600 x 10 seconds, respectively. In hours the times 
are 1.4, 3.3, 22.2 and 167 hours. Problems of this size are 
possible on all but the last machine and can in fact be performed 
on a personal computer (486) if one day turnarounds are 
permissible. 

V. VERSION C CODE CAPABILITIES 

The Penn State University FDTD Electromagnetic Scattering 
Code Version C has the following capabilities: 

1) Ability to model lossy dielectric and perfectly conducting 
scatterers . 

2) Ability to model lossy magnetic scatterers. 

3) First and second order outer radiation boundary condition 
(ORBC) operating on the electric fields for dielectric 
scatterers. 

4) First and second order ORBC operating on the magnetic fields 
for magnetic scatterers. 

5) Near to far zone transformation capability to obtain far zone 
scattered fields. 

6) Gaussian and smooth cosine incident waveforms with arbitrary 
incidence angles. 

7) Near zone field, current or power sampling capability. 

8) Companion code for computing Radar Cross Section (RCS) . 

VI. DEFAULT SCATTERING GEOMETRY 

The code as delivered is set up to calculate the far zone 
backscatter fields for a 20 )l cm radius, lossy magnetic sphere with 
parameters e=e 0 , m=4m 0 and a=2.8413E+4. The problem space size 
is 66 by 66 by 66 cells in the x, y and z directions, the cells 
are 0.01 cm cubes, and the incident waveform is a 0-polarized 
Gaussian pulse with incidence angles of 0=22.5 and 0=22.5 
degrees. The output data files are included as a reference along 
with a code (RCS3D.F0R) for computing the frequency domain RCS 
using these output data files. The ORBC is the second order 
absorbing boundary condition set forth by Mur [2]. 
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VII. SUBROUTINE DESCRIPTION 

In the description for each subroutine, an asterisk (*) will 
be placed by the subroutine name if that particular subroutine is 
normally modified when defining a scattering problem. 

MAIN ROUTINE 

The main routine in the program contains the calls for all 
necessary subroutines to initialize the problem space and 
scattering object (s) and for the incident waveform, far zone 
transformation, field update subroutines, outer radiation 
boundary conditions and field sampling. 

The main routine begins with the include statement and then 
appropriate data files are opened, and subroutines ZERO, BUILD 
and SETUP are called to initialize variables and/or arrays, build 
the object (s) and initialize the incident waveform and 
miscellaneous parameters, respectively. Subroutine SETFZ is 
called to intialize parameters for the near to far zone 
transformation if far zone fields are desired. 

The main loop is entered next, where all of the primary 
field computations and data saving takes place. During each time 
step cycle, the EXSFLD, EYSFLD, and EZSFLD subroutines are called 
to update the x, y, and z components of the scattered electric 
field. The six electric field outer radiation boundary 
conditions (RADE??) are called next to absorb any outgoing 
scattered fields for perfectly conducting or lossy dielectric 
scatterers. Time is then advanced 1/2 time step according to the 
Yee algorithm and then the HXSFLD, HYSFLD, AND HZSFLD subroutines 
are called to update the x, y, and z components of scattered 
magnetic field. The six magnetic field outer radiation boundary 
conditions (RADH??) are called next to absorb any outgoing 
scattered fields for lossy magnetic scatterers. Time is then 
advanced another 1/2 step and then either near zone fields are 
sampled and written to disk in DATSAV, and/or the near zone to 
far zone vector potentials are updated in SAVFZ. The parameter 
NZFZ (described later) in the common file defines the type of 
output fields desired. 

After execution of all time steps in the main field update 
loop, subroutine FAROUT is called if far zone fields are desired 
to compute the far zone fields and write them to disk. At this 
point, the execution is complete. 

SUBROUTINE SETFZ 

This subroutine initializes the necessary parameters 
required for far zone field computations. The code as furnished 
computes backscatter far zone fields and can compute bistatic far 
zone fields for one scattering angle (i.e. one 0 and <p angle) . 
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Refer to reference [3] for a complete description of the near to 
far zone transformation. Other versions of this subroutine 
provide for multiple bistatic angles. 

SUBROUTINE SAVFZ 

This subroutine updates the near zone to far zone vector 
potentials. 

SUBROUTINE FAROUT 

This subroutine changes the near zone to far zone vector 
potentials to far zone electric field 0 and 0 components and 
writes them to disk. 

SUBROUTINE BUILD * 

This subroutine "builds" the scattering object (s) by 
initializing the I DONE , IDTWO, IDTHRE, IDFOR, IDFIV and IDSIX 
arrays. The IDONE-IDTHRE arrays are for specifying perfectly 
conducting and lossy dielectric materials. The IDFOR-IDSIX 
arrays are for lossy magnetic materials. The reason for the 
separate arrays is so the user can independently control the 
exact placement of dielectric and magnetic material in the Yee 
cells. Refer to Figure 1 for a diagram of the basic Yee cell. 

For example, setting an element of the IDONE array at some I,J,K 
location is actually locating dielectric material at a cell edge 
whose center location is I+0.5,J,K. Setting an element of the 
IDFOR array at some I,J,K location is actually locating magnetic 
material perpendicular to a cell face whose center location is 
I , J+0 . 5 , K+0 . 5 , or equivalently, locating magnetic material at an 
edge on the dual H field mesh. The difference between the IDONE 
and IDFOR array locations is a direct result of the field offsets 
in the Yee cell (see Figure 1) . Thus, materials with diagonal 
permittivity and/or diagonal permeability tensors can be modeled. 
The default material type for all ID??? arrays is 0, or free 
space. By initializing these arrays to values other than 0, the 
user is defining an object by determining what material types are 
present at each spatial location. Other material types available 
for IDONE-IDTHRE are 1 for perfectly conducting objects, and 2-9 
for lossy dielectrics. IDONE-IDTHRE are normally set to 0 for 
magnetic scatterers. Other material types available for IDFOR- 
IDSIX are 10-19 for lossy magnetic materials. IDFOR-IDSIX are 
normally set to 0 for perfectly conducting or dielectric 
scatterers. If the user wants a material with both dielectric 
and magnetic properties (i.e. permittivity other than e 0 for 
magnetic materials, and permeability other than /u 0 for dielectric 
materials) , then he/she must define IDONE-IDSIX for that 
particular material. This subroutine also has a section that 
checks the ID??? arrays to determine if legal material types have 
been defined throughout th£ problem space. The actual material 
parameters (e, n, a, and a ) are defined in subroutine SETUP. 
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The default geometry is a 20 cm radius, lossy magnetic sphere. 

The user must be careful that his/her object created in the 
BUILD subroutine is properly formed. There is not a direct 
one-to-one correspondence between the dielectric and magnetic 
ID??? arrays. However, one can define a correspondence, so that 
code used to generate a dielectric object can be modified to 
generate a magnetic object. 

To see this consider that we have set the permittivity at 
cell locations corresponding to 


EX ( I , J , K) , EY ( I , J , K) , EZ ( I , J , K) 

using the IDONE, IDTWO, and IDTHRE arrays respectively. This 
determines one corner of a dielectric cube. If we wish to define 
the corner of a corresponding magnetic cube, offset 1/2 cell in 
the x, y, z directions, we would set the locations of the 
magnetic fields 


HX ( 1+1 , J , K) , HY ( I , J+l , K) , HZ ( I , J, K+l) 

as magnetic material using the IDFOR, IDFIV, and IDSIX arrays. 

This example indicates the following correspondence between 
BUILDing dielectric and magnetic objects: 


Dielectric 

Object 


IDONE (I, J,K) 
IDTWO (I, J,K) 
IDTHRE ( I, J,K) 


Magnetic 

Object 


= IDFOR(I+l, J,K) 
= IDFIV (I, J+l, K) 
= IDSIX ( I , J , K+l) 


To illustrate this for a somewhat more complicated case, consider 
an example of a 3x3 cell dielectric plate located in the XY plane 
at K=5. The plate can be generated with the following FORTRAN 
lines : 

PLATE 


DO 2 0 J=4 , 7 
DO 10 1=4,7 

IF ( I . NE . 7 ) IDONE ( I , J , K) =2 
IF ( J . NE . 7 ) IDTWO ( I , J , K) =2 
10 CONTINUE 
20 CONTINUE 


7 


6 

J 

A 5 


4 1 

4 5 6 7 

>X 
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If the same FORTRAN lines were used to try to generate a magnetic 
plate, the object generated would actually be unconnected: 


7. 


DO 20 J=4 , 7 6 

DO 10 1=4,7 

IF ( I . NE . 7 ) IDFOR (I,J,K)=11 J 

IF ( J. NE. 7 ) IDFIV (I,J,K)=11 - 5 

10 CONTINUE ] 

20 CONTINUE j 

! 4 ' 

! 4 


The correct way to build the magnetic plate is with the following 
FORTRAN code: 

DO 20 J=4 , 7 
DO 10 1=4,7 

IF ( I . NE . 4 ) IDFOR ( I , J, K) =11 
IF ( J. NE . 4 ) IDFIV ( I , J, K) =11 
10 CONTINUE 
20 CONTINUE 

This is equivalent to the correspondence relationship given 
above. See comments in the BUILD subroutine for further 
explanation of the ID??? arrays. 

When it is important to place the object in the center of 
the problem space (to have lowest possible cross-pol scattering 
for symmetric objects) , NX etc. should be odd for dielectrics and 
even for magnetics. This is due to the field locations in the 
Yee cell and also the placement of the E and H field absorbing 
boundary condition surfaces. 

If the object being modeled has curved surfaces, edges, etc. 
that are at an angle to one or more of the coordinate axes, then 
that shape must be approximately modeled by lines and faces in a 
"stair-stepped" (or stair-cased) fashion. This stair-cased 
approximation introduces errors into computations at higher 
frequencies. Intuitively, the error becomes smaller as more 
cells are used to stair-case a particular object. Thus, the 
default Nickel Ferrite sphere scattering geometry is a stair- 
cased sphere . 

When the user's test object is dielectric it is apparently 
better to use the Mur E field boundary formulation. For a 
magnetic scattering object it seems better to use the Mur H field 
boundary formulation. The reason for this has to do with the 
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reactive part of the energy (near zone fields) reaching the 
absorbing boundary surface. These near zone fields differ for 
magnetic and dielectric objects, even for dual objects where the 
far fields are essentially identical. 

One final comment on building dual dielectric vs magnetic 
objects can be explained by considering the duality between 
magnetic and dielectric materials. In testing the code, it was 
helpful to test dual dielectric/magnetic objects since they 
scatter identically in the far zone. In specifying the material 
of a magnetic scatterer to be the dual of a dielectric, the 
following transformation must be applied: 

DIELECTRIC < — DUAL — > MAGNETIC 


E field < 

H field < 


e < 

M < 

6 (k) < 

a < 


> H field 

> -E field 

> M 

> e 

> B (k) 

> tf*(M 0 / € o) =a 


The somewhat surprising entry in this table is the relationship 
between the dielectric and magnetic conductivities. The reason 
why the magnetic conductivity, a , for magnetic materials has to 
be scaled as indicated is that for duality to be applied the free 
space impedance must be inverted. Although this is not generally 
possible for actual problems, identical far zone scattering for 
dielectric apd magnetic scatterers can be achieved by this 
scaling of a as above and realizing that E and H scattered £ield 
to incident field ratios will not invert. This scaling of a is 
why conductivity for magnetic materials is usually not a 
dominating feature, and in fact is often neglected. 


SUBROUTINE DCUBE 


This subroutine builds cubes of dielectric material by 
defining four each of I DONE , IDTWO and IDTHRE components 
corresponding to one spatial cube of dielectric material. It can 
also be used to define thin (i.e. up to one cell thick) 
dielectric or perfectly conducting plates. Refer to comments 
within DCUBE for a description of the arguments and usage of the 
subroutine. This subroutine could be modified to build cubes 
and/or plates of magnetic materials by using a triple do-loop in 
BUILD (after calls to DCUBE) over coordinate indices I, J, K and 
applying the correspondence between the IDONE-IDTHRE and IDFOR- 
IDSIX arrays that was previously discussed. 

SUBROUTINE SETUP * 


This subroutine initializes many of the constants required 
for incident field definition, field update equations, outer 
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radiation boundary conditions and # material parameters. The 
material parameters e, /x, a and a are defined for each material 
type using the material arrays EPS, XMU, SIGMA and SIGMAC 
respectively. The array EPS is used for the total permittivity, 
XMU is used for the total permeability, SIGMA is used for the 
electric conductivity and SIGMAC is used for the magnetic 
conductivity (useful for running dual problems) . These arrays 
are initialized in SETUP to free space material parameters for 
all material types and then the user is required to modify these 
arrays for his/her scattering materials. Thus, for the lossy 
dielectric material type 2, the user must define EPS (2) and 
SIGMA (2). The remainder of the subroutine computes constants 
used in field update equations and boundary conditions and writes 
the diagnostics file. 

SUBROUTINE EXSFLD 

This subroutine updates all x components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDONE array are used to determine the type of material present 
and the corresponding update equation to be used. These 
scattered field equations are based on the development given in 
[4]. 

SUBROUTINE EYSFLD 

This subroutine updates all y components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDTWO array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINE EZSFLD 

This subroutine updates all z components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDTHRE array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINES RADEYX, RADEZX, RADEZY, RADEXY, RAD EX Z and RADEYZ 

These subroutines apply the outer radiation boundary 
conditions to the scattered electric field on the outer 
boundaries of the problem space for non-magnetic scatterers. The 
user controls selection of these routines through the parameter 
MAGNET (described later) . 
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SUBROUTINE HXSFLD 

This subroutine updates all x components of scattered 
magnetic field at each time step. IF statements based upon the 

IDFOR array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINE HYSFLD 

This subroutine updates all y components of scattered 
magnetic field at each time step. IF statements based upon the 
IDFIV array are used to determine the type of material present 
and the corresponding update eguation to be used. 

SUBROUTINE HZSFLD 

This subroutine updates all z components of scattered 
magnetic field at each time step. IF statements based upon the 
IDSIX array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINES RADHXZ , RADHYX, RADHZY, RADHXY , RADHYZ and RADHZX 

These subroutines apply the outer radiation boundary 
conditions to the scattered magnetic field on the outer 
boundaries of the problem space for magnetic scatterers. The 
user controls selection of these routines through the parameter 
MAGNET (described later) . 

SUBROUTINE DATSAV * 

This subroutine samples near zone scattered field quantities 
and saves them to disk. This subroutine is where the quantities 
to be sampled and their spatial locations are to be specified and 
is only called if near zone fields only are desired or if both 
near and far zone fields are desired. Total field quantities can 
also be sampled. See comments within the subroutine for 
specifying sampled scattered and/or total field quantities. When 
sampling magnetic fields, remember the St/2 time difference 
between E and H when writing the fields to disk. Sections of 
code within this subroutine determine if the sampled quantities 
and the spatial locations have been properly defined. 

FUNCTIONS EXI, EYI , EZI, HXI , HYI and HZI 

These functions are called to compute the x, y and z 
components of incident electric and magnetic field, respectively. 
The functional form of the incident field is contained in a 
separate function SOURCE. 



16 


FUNCTION SOURCE * 

This function contains the functional form of the incident 
field. The code as furnished uses the Gaussian form of the 
incident field. An incident smooth cosine pulse is also 
available by uncommenting the required lines and commenting out 
the Gaussian pulse. Thus, this function need only be modified if 
the user changes the incident pulse from Gaussian to smooth 
cosine. A slight improvement in computing speed and 
vectorization may be achieved by moving this function inside each 
of the incident field functions EXI, EYI and so on. 

FUNCTIONS DEXI, DEYI , DEZI, DHXI , DHYI and DHZI 

These functions are called to compute the x, y and z 
components of the time derivative of incident electric and 
magnetic field, respectively. The functional form of the 
incident field is contained in a separate function DSRCE. 

FUNCTION DSRCE * 

This function contains the functional form of the time 
derivative of the incident field. The code as furnished uses the 
time derivative of the Gaussian form of the incident field. A 
smooth cosine pulse time derivative is also available by 
uncommenting the required lines and commenting out the Gaussian 
pulse. Thus, the function need only be modified if the user 
changes from the Gaussian to smooth cosine pulse. Again, a 
slight improvement in computing speed and vectorization may be 
achieved by moving this function inside each of the time 
derivative incident field functions DEXI, DEYI and so on. 

SUBROUTINE ZERO 

This subroutine initializes various arrays and variables to 

zero. 

VIII. INCLUDE FILE DESCRIPTION (COMMONC.FOR) * 

The include file, COMMONC.FOR, contains all of the arrays 
and variables that are shared among the different subroutines. 
This file will require the most modifications when defining 
scattering problems. A description of the parameters that are 
normally modified follows. 

The parameters NX, NY and NZ specify the size of the problem 
space in cells in the x, y and z directions respectively. For 
problems where it is crucial to center the object within the 
problem space, then NX, NY and NZ should be odd for dielectric 
scatterers and even for magnetic scatterers. The parameter NTEST 
defines the number of near zone quantities to be sampled and NZFZ 
defines the field output format. Set NZFZ=0 for near zone fields 
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only, NZFZ=1 for far zone fields only and NZFZ=2 for both near 
and far zone fields. Parameter MAGNET is used to define magnetic 
scatterers and it controls the choice of RADE?? versus RADH?? 
absorbing boundary subroutines. It is set to 0 for dielectric 
scatterers and to 1 for magnetic scatterers. Parameter NSTOP 
defines the maximum number of time steps. DELX, DELY, and DELZ 
(in meters) define the cell size in the x, y and z directions 
respectively. The 0 and 0 incidence angles (in degrees) are 
defined by THINC and PHINC respectively and the polarization is 
defined by ETHINC and EPHINC. ETHINC=1.0, EPHINC=0.0 for 0 - 
polarized incident field and ETHINC=0.0, EPHINC=1.0 for 0- 
polarized incident fields. Parameters AMP and BETA define the 
maximum amplitude and the temporal width of the incident 
pulse respectively. BETA automatically adjusts when the cell 
size is changed and normally should not be changed by the user. 
The far zone scattering angles are defined by THETFZ and PHIFZ. 
The code as furnished performs backscatter computations, but 
these parameters could be modified for a bistatic computation. 

IX. RCS COMPUTATIONS 

A companion code, RCS3D.F0R, has been included to compute 
RCS versus frequency. It uses the file name of the FDTD far zone 
output data (FZ0UT3D.DAT) and writes data files of far zone 
electric fields versus time (FZTIME.DAT) and RCS versus frequency 
(3DRCS.DAT). The RCS computations are performed up to the 10 
cell/l 0 frequency limit. Refer to comments within this code for 
further details. 

X. RESULTS 

As previously mentioned, the code as furnished models a 20 
cm radius, lossy magnetic sphere and computes backscatter far 
zone scattered fields at angles of 0=22.5 and 0=22.5 degrees. 
Results are included for the dual dielectric sphere and the 
default magnetic sphere. The material parameters for the 
dielectric dual of the magnetic material are e=4e 0 , a=0.2 and 
/x=/x 0 . For these materials there are 5 cells per wavelength at 
approximately 3.0 GHz. 

Figures 2-3 show the co-polarized far zone scattered field 
versus time and the co-polarized RCS versus frequency for the 
dielectric sphere. 

Figures 4-5 show the co-polarized far zone electric field 
versus time and the co-polarized RCS for the magnetic sphere. 

XI. SAMPLE PROBLEM SETUP 

The code as furnished models a 20 cm radius, lossy magnetic 
sphere and computes backscatter far zone scattered fields at 
angles of 0=22.5 and 0=22.5 degrees. The corresponding output 
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data files are also provided, along with a code to compute Radar 
Cross Section using these data files. In order to change the 
code to a new problem, many different parameters need to be 
modified. A sample problem setup will now be discussed. 

Suppose that the problem to be studied is RCS backscatter 
versus frequency from a 28 cm by 31 cm perfectly conducting plate 
with a 3 cm dielectric coating with a dielectric constant of 4e 0 
using a 0-polarized field. The backscatter angles are 0=30.0 and 
0=60.0 degrees and the frequency range is up to 3 Ghz. 

Since the frequency range is up to 3 Ghz, the cell size must 
be chosen appropriately to resolve the field IN ANY MATERIAL at 
the highest frequency of interest. A general rule is that the 
cell size should be 1/10 of the wavelength at the highest 
frequency of interest. For difficult geometries, 1/20 of a 
wavelength may be necessary. The free space wavelength at 3 GHz 
is A 0 =10 cm and the wavelength in the dielectric coating at 3 GHz 
is 5 cm. The cell size is chosen as 1 cm, which provides a 
resolution of 5 cells/A in the dielectric coating and 10 cells/A 0 
in free space. Numerical studies have shown that choosing the 

cell size < 1/4 of the shortest wavelength in any material is 
the practical lower limit. Thus the cell size of 1 cm is barely 
adequate. The cell size in the x, y and z directions is set in 
the common file through variables DELX, DELY and DELZ. Next the 
problem space size must be large enough to accomodate the 
scattering object, plus at least a five cell boundary (10 cells 
is more appropriate) on every side of the object to allow for the 
far zone field integration surface. It is advisable for plate 
scattering to have the plate centered in the x and y directions 
of the problem space in order to reduce the cross-polarized 
backscatter and to position the plate low in the z direction to 
allow strong specular reflections multiple encounters with the 
ORBC . A 10 cell border is chosen, and the problem space size is 
chosen as 49 by 52 by 49 cells in the x, y and z directions 
respectively. As an initial estimate, allow 2048 time steps so 
that energy trapped within the dielectric layer will radiate. 

Thus parameters NX, NY and NZ in COMMONC.FOR would be changed to 
reflect the new problem space size, and parameter NSTOP is 
changed to 2048. If all transients have not been dissipated 
after 2048 time steps, then NSTOP will have to be increased. 
Truncating the time record before all transients have dissipated 
will corrupt frequency domain results. Parameter NZFZ must be 
equal to 1 since we are interested in far zone fields only. 
Parameter MAGNET must be equal to 0 for the dielectric scatterer. 
To build the object, the following lines are inserted into the 
BUILD subroutine: 

C 

C BUILD THE DIELECTRIC SLAB FIRST 

C 


ISTART=1 1 
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JSTART=11 

KSTART=11 

NXWIDE=28 

NYWIDE=3 1 

NZWIDE=3 

MTYPE=2 

CALL DCUBE (ISTART , JSTART , K START , NXWIDE , NYWIDE , NZWIDE , MTYPE ) 
C 

C BUILD PEC PLATE NEXT 

C 

ISTART=11 

JSTART=11 

KSTART=11 

NXWIDE=28 

NYWIDE=3 1 

NZWIDE=0 

MTYPE=1 

CALL DCUBE (ISTART , JSTART , KSTART, NXWIDE , NYWIDE , NZWIDE , MTYPE) 

The PEC plate is built last on the bottom of the dielectric 
slab to avoid any air gaps between the dielectric material and 
the PEC plate. In the common file, the incidence angles THINC 
and PHINC have to be changed to 30.0 and 60.0 respectively, the 
cell sizes (DELX, DELY, DELZ) are set to 0.01, and the 
polarization is set to ETHINC=1.0 and EPHINC=0.0 for 0-polarized 
fields. Since dielectric material 2 is being used for the 
dielectric coating, the constitutive parameters EPS (2) and 
SIGMA (2) are set to 4e 0 and 0.0 respectively, in subroutine 
SETUP. This completes the code modifications for the sample 
problem. 

XII. NEW PROBLEM CHECKLIST 

This checklist provides a quick reference to determine if 
all parameters have been defined properly for a given scattering 
problem. A reminder when defining quantities within the code: 
use MKS units and specify all angles in degrees. 

COMMONC . FOR : 

1) Is the problem space sized correctly? (NX, NY, NZ) 

2) For near zone fields, is the number of sample points correct? 
(NTEST) 

3) Is parameter NZFZ defined correctly for desired field 
outputs? 

4) Is parameter MAGNET defined correctly for the type of 
scatterer? 

5) Is the number of time steps correct? (NSTOP) 
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6) Are the cell dimensions (DELX, DELY, DELZ) defined correctly? 

7) Are the incidence angles (THINC, PHINC) defined correctly? 

8) Is the polarization of the incident wave defined correctly 
(ETHINC, EPHINC)? 

9) For other than backscatter far zone field computations, are 
the scattering angles set correctly? (THETFZ, PHIFZ) 

SUBROUTINE BUILD: 

1) Is the object completely and correctly specified? 

SUBROUTINE SETUP: 

1) Are the constitutive parameters for each material specified 
correctly? (EPS, XMU, SIGMA, SIGMAC) 

FUNCTIONS SOURCE and DSRCE: 

1) If the Gaussian pulse is not desired, is it commented out and 
the smooth cosine pulse uncommented? 

SUBROUTINE DATSAV : 

1) For near zone fields, are the sampled field types and spatial 
locations correct for each sampling point? (NTYPE, IOBS, JOBS, 
KOBS) 

XIII. REFERENCES 

[1] K. S. Yee, "Numerical solution of initial boundary value 
problems involving Maxwell's equations in isotropic media," 
IEEE Trans. Antennas Propaaat. . vol. AP-14, pp. 302-307, May 
1966. 

[2] G. Mur, "Absorbing boundary conditions for the Finite- 
Difference approximation of the Time-Domain Electromagnetic- 
Field Equations," IEEE Trans. Electromaan. Compat. . vol. 

EMC- 2 3 , pp. 377-382, November 1981. 

[3] R. J. Luebbers et. al., "A Finite Difference Time-Domain 
Near Zone to Far Zone Transformation," IEEE Trans. Antennas 
Propaqat. . vol. AP-39, no. 4, pp. 429-433, April 1991. 

[4] R. Holland, L. Simpson and K. S. Kunz, "Finite-Difference 
Time-Domain Analysis of EMP Coupling to Lossy Dielectric 
Structures," IEEE Trans. Electromaan. Compat. . vol. EMC-22, 
pp. 203-209, August 1980. 



XIV. 

Fig. 

Fig. 

Fig. 

Fig. 

Fig. 


21 


FIGURE TITLES 

1 Standard three dimensional Yee cell showing placement 
of electric and magnetic fields. 

2 Co-polarized far zone scattered field versus time for 
dielectric sphere with 20 cm radius. 

3 Co-polarized Radar Cross Section versus frequency for 
dielectric sphere with 20 cm radius. 

4 Co-polarized far zone scattered field versus time for 
magnetic sphere with 20 cm radius. 

5 Co-polarized Radar Cross Section versus frequency for 
magnetic sphere with 20 cm radius. 
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I . INTRODUCTION 

The Penn State Finite Difference Time Domain Electromagnetic 
Scattering Code Version D is a three dimensional numerical 
electromagnetic scattering code based upon the Finite Difference 
Time Domain Technigue (FDTD) . The supplied version of the code 
is one version of our current three dimensional FDTD code set. 
This manual provides a description of the code and corresponding 
results for several scattering problems. The manual is organized 
into fourteen sections: introduction, description of the FDTD 

method, operation, resource requirements. Version D code 
capabilities, a brief description of the default scattering 
geometry, a brief description of each subroutine, a description 
of the include file (COMMOND. FOR) , a section briefly discussing 
Radar Cross Section (RCS) computations, a section discussing some 
scattering results, a sample problem setup section, a new problem 
checklist, references and figure titles. 

II. FDTD METHOD 

The Finite Difference Time Domain (FDTD) technique models 
transient electromagnetic scattering and interactions with 
objects of arbitrary shape and/or material composition. The 
technique was first proposed by Yee [1] for isotropic, non- 
dispersive materials in 1966; and has matured within the past 
twenty years into a robust and efficient computational method. 

The present FDTD technique is capabable of transient 
electromagnetic interactions with objects of arbitrary and 
complicated geometrical shape and material composition over a 
large band of frequencies. This technique has recently been 
extended to include dispersive dielectric materials, chiral 
materials and plasmas. 

In the FDTD method, Maxwell's curl equations are discretized 
in time and space and all derivatives (temporal and spatial) are 
approximated by central differences. The electric and magnetic 
fields are interleaved in space and time and are updated in a 
second-order accurate leapfrog scheme. The computational space 
is divided into cells with the electric fields located on the 
edges and the magnetic fields on the faces (see Figure 1) . FDTD 
objects are defined by specifying dielectric and/or magnetic 
material parameters at electric and/or magnetic field locations. 

Two basic implementations of the FDTD method are widely used 
for electromagnetic analysis: total field formalism and scattered 
field formalism. In the total field formalism, the electric and 
magnetic field are updated based upon the material type present 
at each spatial location. In the scattered field formalism, the 
incident waveform is defined analytically and the scattered field 
is coupled to the incident field through the different material 
types. For the incident field, any waveform, angle of incidence 
and polarization is possible. The separation of the incident and 
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scattered fields conveniently allows an absorbing boundary to be 
employed at the extremities of the discretized problem space to 
absorb the scattered fields. 

This code is a scattered field code, and the total E and H 
fields may be found by combining the incident and scattered 
fields. Any type of field quantity (incident, scattered, or 
total) , Poynting vector or current are available anywhere within 
the computational space. These fields, incident, scattered and 
total, may be found within, on or about the interaction object 
placed in the problem space. By using a near to far field 
transformation, far fields can be determined from the near fields 
within the problem space thereby affording radiation patterns and 
RCS values. The accuracy of these calculations is typically 
within a dB of analytic solutions for dielectric and magnetic 
sphere scattering. Further improvements are expected as better 
absorbing boundary conditions are developed and incorporated. 

III. OPERATION 


Typically, a truncated Gaussian incident waveform is used to 
excite the system being modeled, however certain code versions 
also provide a smooth cosine waveform for convenience in modeling 
dispersive materials. The interaction object is defined in the 
discretized problem space with arrays at each cell location 
created by the discretization. All three dielectric material 
types for E field components within a cell can be individually 
specified by the arrays IDONE(I, J,K) , IDTWO(I, J,K) , 

IDTHRE(I, J,K) . This models arbitrary dielectric materials with n 
= By an obvious extension to six arrays, magnetic materials 

with fi ^ Hq can be modeled. 

Scattering occurs when the incident wave, marched forward in 
time in small steps set by the Courant stability condition, 
reaches the interaction object. Here a scattered wave must 
appear along with the incident wave so that the Maxwell equations 
are satisfied. If the material is a perfectly conductive metal 
then only the well known boundary condition 


must be satisfied. For a nondispersive dielectric the 
requirement is that the total field must satisfy the Maxwell 
equations in the material: 


Vx E tot 


Vx ( E inc + E scat 


1 dH tot 
M 0 dt 


1 3(H inc + H scat ) 
/i dt 

'o 


(2) 
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VxH tot = Vx (H ,nc +H scat ) 


3E tot 
= e 

at 


+ aE 


tot 


(3) 


a ( e' oc + E scat ) 

= e 

at 


+ <7 (E inc + E scat ) 


(4) 


Additionally the incident wave, defined as moving unimpeded 
through a vacuum in the problem space, satisfies everywhere in 
the problem the Maxwell equations for free space 


Vx E inc 


1 aH inc 
dt 


(5) 


VxH ,nc = e 


aE' 


at 


( 6 ) 


Subtracting the second set of equations from the first yields the 
Maxwell equations governing the scattered fields in the material: 


Vx E scat 


i aH scat 
dt 


(7) 


a E ,nc 

VxH scat =(£-e ) +aE ,nc +e 

0 at at 


^T? SCat 

atL ,~scat 

+aE 


( 8 ) 


Outside the material this simplifies to: 


Vx E 


scat 


1 3H 


scat 


M 0 dt 


(9) 


VxH scat = c 


aE : 


scat 


at 


( 10 ) 


Magnetic materials, dispersive effects, non-linearities, 
etc. , are further generalizations of the above approach. Based 
on the value of the material type, the subroutines for 
calculating scattered E and H field components branch to the 
appropriate expression for that scattered field component and 
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that component is advanced in time according to the selected 
algorithm. As many materials can be modeled as desired, the 
number equals the dimension selected for the flags. If materials 
with behavior different from those described above must be 
modeled, then after the appropriate algorithm is found, the 
code's branching structure allows easy incorporation of the new 
behavior. 

IV. RESOURCE REQUIREMENTS 

The number of cells the problem space is divided into 
the six components per cell set the problem space storage 
requirements 

Storage = NC x 6 components /cell x 4 bytes /component 

and the computational cost 

Operations = NC x 6 comp/cell x 10 ops/component x N (12) 

where N is the number of time steps desired. 

N typically is on the order of ten times the number of cells 
on one side of the problem space. More precisely for cubical 

cells it takes time steps to traverse a single cell when the 
time step is set by the Courant stability condition 

Ax 

At = Ax = cell size dimension ( 13 ) 

ho 


The condition on N is then that 

1 

7 number cells on a side ( 14 ) 

NC ~ 

of the problem space 

The earliest aircraft modeling using FDTD with approximately 30 
cells on a side required approximately 500 time steps. For more 
recent modeling with approximately 100 cells on a side, 2000 or 
more time steps are used. 

For (100 cell) 3 problem spaces, 24 MBytes of memory are 
required to store the fields. Problems on the order of this size 
have been run on a Silicon Graphics 4D 220 with 32 MBytes of 
memory, IBM RISC 6000, an Intel 486 based machine, and VAX 
11/785. Storage is only a problem as in the case of the 486 
where only 16 MBytes of memory was available. 3 This limited the 
problem space size to approximately (80 cells) . 


N 


10 x ( 


/Tn 


C J ) 


times 

an 


G-J. 
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For (100 cell) 3 problems with approximately 2000 time steps, 
there is a total of 120 x 10 7 8 9 operations to perform. The speeds 
of the previously mentioned machines are 24 MFLOPs (4 processor 
upgraded version), 10 MFLOPS, 1.5 MFLOPS, ^nd 0.2 MFLOPs. T^e 
run times are then 5 x 10 3 seconds, 12 x 10 seconds, 80 x 10 
seconds and 600 x 10 3 seconds, respectively. In hours the times 
are 1.4, 3.3, 22.2 and 167 hours. Problems of this size are 
possible on all but the last machine and can in fact be performed 
on a personal computer (486) if one day turnarounds are 
permissible. 

V. VERSION D CODE CAPABILITIES 

The Penn State University FDTD Electromagnetic Scattering 
Code Version D has the following capabilities: 

1) Ability to model lossy dielectric and perfectly conducting 
scatterers . 

2) Ability to model lossy magnetic scatterers. 

3) Ability to model dispersive dielectric and dispersive 
magnetic scatterers. This dispersive FDTD method is now 
designated (FD)TD for Frequency-Dependent Finite Difference Time 
Domain. 

4) First and second order outer radiation boundary condition 
(ORBC) operating on the electric fields for dielectric and 
dispersive dielectric scatterers. 

5) First and second order ORBC operating on the magnetic fields 
for magnetic and dispersive magnetic scatterers. 

6) Near to far zone transformation capability to obtain far zone 
scattered fields. 

7) Gaussian and smooth cosine incident waveforms with arbitrary 
incidence angles. 

8) Near zone field, current or power sampling capability. 

9) Companion code for computing Radar Cross Section (RCS) . 

VI. DEFAULT SCATTERING GEOMETRY 

The code as delivered is set up to calculate the far zone 
backscatter fields for a 6.67 meter radius, dispersive, Nickel 
Ferrite sphere. Nickel Ferrite is defined by a frequency 
dependent permeability given by 

where i± 0 is the infinite frequency permeability, n s is the zero 



9 


^0 


M s -M. 

= + — 

i + :o>r 0 


(15) 


frequency permeability, r Q is the relaxation time and g> is the 
radian frequency. The Nickel Ferrite parameters are n o =l, /x s =100 
and r Q = 20 ns. The problem space size is 66 by 66 by 66 cells in 

the x, y and z directions, the cells are 1/3 m cubes, and the 
incident waveform is a 0— polarized smooth cosine pulse with 
incidence anqles of 0=22.5 and 0=22.5 degrees. The output data 
files are included as a reference along with a code (RCS3D.FOR) 
for computing the frequency domain RCS using these output data 
files. The ORBC is the second order absorbing boundary condition 
set forth by Mur [ 2 ] . 

VII. SUBROUTINE DESCRIPTION 

In the description for each subroutine, an asterisk (*) will 
be placed by the subroutine name if that particular subroutine is 
normally modified when defining a scattering problem. 

MAIN ROUTINE 

The main routine in the program contains the calls for all 
necessary subroutines to initialize the problem space and 
scattering object (s) and for the incident waveform, far zone 
transformation, field update subroutines, outer radiation 
boundary conditions and field sampling. 

The main routine begins with the include statement and then 
appropriate data files are opened, and subroutines ZERO, BUILD 
and SETUP are called to initialize variables and/or arrays, build 
the object (s) and initialize the incident waveform and 
miscellaneous parameters, respectively. Subroutine SETFZ is 
called to intialize parameters for the near to far zone 
transformation if far zone fields are desired. 

The main loop is entered next, where all of the primary 
field computations and data saving takes place. During each time 
step cycle, the EXSFLD , EYSFLD, and EZSFLD subroutines are called 
to update the x, y, and z components of the scattered electric 
field. The six electric field outer radiation boundary 
conditions (RADE??) are called next to absorb any outgoing 
scattered fields for perfectly conducting, dielectric, or 
dispersive dielectric scatterers. Time is then advanced 1/2 time 
step according to the Yee algorithm and then the HXSFLD, HYSFLD, 
AND HZSFLD subroutines are called to update the x, y, and z 
components of scattered magnetic field. The six magnetic field 
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outer radiation boundary conditions (RADH??) are called next to 
absorb any outgoing scattered fields for magnetic or dispersive 
magnetic scatterers. Time is then advanced another 1/2 step and 
then either near zone fields are sampled and written to disk in 
DATSAV, and/or the near zone to far zone vector potentials are 
updated in SAVFZ. The parameter NZFZ (described later) in the 
common file defines the type of output fields desired. 

After execution of all time steps in the main field update 
loop, subroutine FAROUT is called if far zone fields are desired 
to compute the far zone fields and write them to disk. At this 
point, the execution is complete. 

SUBROUTINE SETFZ 

This subroutine initializes the necessary parameters 
required for far zone field computations. The code as furnished 
computes backscatter far zone fields and can compute bistatic far 
zone fields for one scattering angle (i.e. one 0 and <p angle) . 
Refer to reference [3] for a complete description of the near to 
far zone transformation. Other versions of this subroutine 
provide for multiple bistatic angles. 

SUBROUTINE SAVFZ 

This subroutine updates the near zone to far zone vector 
potentials. 

SUBROUTINE FAROUT 

This subroutine changes the near zone to far zone vector 
potentials to far zone electric field 0 and <p components and 
writes them to disk. 

SUBROUTINE BUILD * 

This subroutine "builds" the scattering object (s) by 
initializing the IDONE, IDTWO, IDTHRE, IDFOR, IDFIV and IDSIX 
arrays. The I DONE- IDTHRE arrays are for specifying perfectly 
conducting, lossy dielectrics and dispersive dielectric 
materials. The IDFOR-IDSIX arrays are for lossy magnetic and 
dispersive magnetic materials. The reason for the separate 
arrays is so the user can independently control the exact 
placement of dielectric and magnetic material in the Yee cells. 
Refer to Figure 1 for a diagram of the basic Yee cell. For 
example, setting an element of the IDONE array at some I,J,K 
location is actually locating dielectric material at a cell edge 
whose center location is I+0.5,J,K. Setting an element of the 
IDFOR array at some I,J,K location is actually locating magnetic 
material perpendicular to a cell face whose center location is 
I , J+0 . 5 , K+o . 5 , or equivalently, locating magnetic material at an 
edge on the dual H field mesh. The difference between the IDONE 
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and IDFOR array locations is a direct result of the field offsets 
in the Yee cell (see Figure 1) . Thus, materials with diagonal 
permittivity and/or diagonal permeability tensors can be modeled. 
The default material type for all ID??? arrays is 0, or free 
space. By initializing these arrays to values other than 0, the 
user is defining an object by determining what material types are 
present at each spatial location. Other material types available 
for IDONE-IDTHRE are 1 for perfectly conducting objects, 2-9 for 
lossy non-magnetic dielectrics, 20-29 for dispersive dielectrics. 
IDONE-IDTHRE are normally set to 0 for magnetic scatterers. Other 
material types available for IDFOR-IDSIX are 10-19 for lossy 
magnetic materials and 30—39 for dispersive magnetic materials. 
IDFOR-IDSIX are normally set to 0 for perfectly conducting or 
dielectric scatterers. If the user wants a material with both 
dielectric and magnetic properties (i.e. permittivity other than 
e 0 for magnetic materials, and permeability other than /x 0 for 
dielectric materials), then he/she must define IDONE-IDSIX for 
that particular material. This subroutine also has a section 
that checks the ID??? arrays to determine if legal material types 
have been defined throughout the problem space,. The actual non- 
dispersive material parameters (e, /i, a , and a ) are defined in 

subroutine SETUP. The dispersive material parameters (e s , e n , 

t q , a, n s , / t q , and cr ) are also defined in a separate section 

in SETUP. The default geometry is a 6.67 m radius, dispersive. 
Nickel Ferrite sphere. 

The user must be careful that his/her object created in the 
BUILD subroutine is properly formed. There is not a direct 
one-to-one correspondence between the dielectric and magnetic 
ID??? arrays. However, one can define a correspondence, so that 
code used to generate a dielectric object can be modified to 
generate a magnetic object. 

To see this consider that we have set the permittivity at 
cell locations corresponding to 

EX ( I , J , K) , EY ( I , J , K) , EZ(I,J,K) 

using the I DONE, IDTWO, and IDTHRE arrays respectively. This 
determines one corner of a dielectric cube. If we wish to define 
the corner of a corresponding magnetic cube, offset 1/2 cell in 
the x, y, z directions, we would set the locations of the 
magnetic fields 

HX ( 1+1 , J , K) , HY ( I , J+l , K) , HZ(I,J,K+1) 
as magnetic material using the IDFOR, IDFIV, and IDSIX arrays. 

This example indicates the following correspondence between 
BUILDing dielectric and magnetic objects: 
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Dielectric Magnetic 

Object Object 

IDONE ( I , J , K) = IDFOR ( 1+1 , J, K) 

IDTWO ( I , J , K) = IDFIV ( I , J+l , K) 

IDTHRE (I , J, K) = IDSIX ( I , J, K+l) 

To illustrate this for a somewhat more complicated case, consider 
an example of a 3x3 cell dielectric plate located in the XY plane 
at K=5. The plate can be generated with the following FORTRAN 
lines: 

7 


6 

J 

A 5 


4 

4 5 6 7 

>1 

If the same FORTRAN lines were used to try to generate a magnetic 
plate, the object generated would actually be unconnected: 

PLATE 

7. • 


DO 20 J=4 , 7 
DO 10 1=4,7 

IF ( I . NE . 7 ) IDONE ( I , J , K) =2 
IF ( J . NE . 7 ) IDTWO ( I , J, K) =2 
10 CONTINUE 
20 CONTINUE 


PLATE 


i 


DO 20 J=4 , 7 
DO 10 1=4,7 

IF (I . NE. 7 ) IDFOR ( I , J, K) =11 
IF ( J. NE . 7 ) IDFIV ( I , J , K) =11 
10 CONTINUE 
20 CONTINUE 




The correct way to build the magnetic plate is with the following 
FORTRAN code: 


DO 20 J=4 , 7 
DO 10 1=4,7 

IF ( I . NE . 4 ) IDFOR ( I , J , K) =11 

IF ( J . NE . 4 ) IDFIV ( I , J , K) =11 
10 CONTINUE 
20 CONTINUE 
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This is equivalent to the correspondence relationship given 
above. See comments in the BUILD subroutine for further 
explanation of the ID??? arrays. 

When it is important to place the object in the center of 
the problem space (to have lowest possible cross-pol scattering 
for symmetric objects) , NX etc. should be odd for dielectrics and 
even for magnetics. This is due to the field locations in the 
Yee cell and also the placement of the E and H field absorbing 
boundary condition surfaces. 

If the object being modeled has curved surfaces, edges, etc. 
that are at an angle to one or more of the coordinate axes, then 
that shape must be approximately modeled by lines and faces in a 
"stair-stepped" (or stair-cased) fashion. This stair-cased 
approximation introduces errors into computations at higher 
frequencies. Intuitively, the error becomes smaller as more 
cells are used to stair-case a particular object. Thus, the 
default Nickel Ferrite sphere scattering geometry is a stair- 
cased sphere. 

When the user’s test object is dielectric it is apparently 
better to use the Mur E field boundary formulation. For a 
magnetic scattering object it seems better to use the Mur H field 
boundary formulation. The reason for this has to do with the 
reactive part of the energy (near zone fields) reaching the 
absorbing boundary surface. These near zone fields differ for 
magnetic and dielectric objects, even for dual objects where the 
far fields are essentially identical. 

One final comment on building dual dielectric vs magnetic 
objects can be explained by considering the duality between 
magnetic and dielectric materials. In testing the code, it was 
helpful to test dual dielectric/magnetic objects since they 
scatter identically in the far zone. In specifying the material 
of a magnetic scatterer to be the dual of a dielectric, the 
following transformation must be applied: 

DIELECTRIC < — DUAL — > MAGNETIC 


E field 
H field 
e 
M 

ft (k) 
a 


< > H field 

< > -E field 

< > M 

< > € 

< > 6 (k) 

< > a . (/i 0 /e 0 )=a 


The somewhat surprising entry in this table is the relationship 
between the dielectric and magnetic conductivities. The reason 
why the magnetic conductivity, a , for magnetic materials has to 
be scaled as indicated is that for duality to be applied the free 
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space impedance must be inverted. Although this is not generally 
possible for actual problems, identical far zone scattering for 
dielectric apd magnetic scatterers can be achieved by this 
scaling of a as above and realizing that E and H scattered £ield 
to incident field ratios will not invert. This scaling of a is 
why conductivity for magnetic materials is usually not a 
dominating feature, and in fact is often neglected. 

SUBROUTINE DCUBE 

This subroutine builds cubes of dielectric material by 
defining four each of I DONE, IDTWO and IDTHRE components 
corresponding to one spatial cube of dielectric material. It can 
also be used to define thin (i.e. up to one cell thick) 
dielectric or perfectly conducting plates. Refer to comments 
within DCUBE for a description of the arguments and usage of the 
subroutine. This subroutine could be modified to build cubes 
and/or plates of magnetic materials by using a triple do-loop in 
BUILD (after calls to DCUBE) over coordinate indices I, J, K and 
applying the correspondence between the IDONE-IDTHRE and IDFOR- 
IDSIX arrays that was previously discussed. 

SUBROUTINE SETUP * 

This subroutine initializes many of the constants required 
for incident field definition, field update equations, outer 
radiation boundary conditions and # material parameters. The 
material parameters e, n, a and a are defined for each material 
type (non-dispersive) using the material arrays EPS, XMU, SIGMA 
and SIGMAC respectively. The array EPS is used for the total 
permittivity, XMU is used for the total permeability, SIGMA is 
used for the electric conductivity and SIGMAC is used for the 
magnetic conductivity (useful for running dual problems) . These 
arrays are initialized in SETUP to free space material parameters 
for all material types and then the user is required to modify 
these arrays for his/her scattering materials. Thus, for the 
lossy dielectric material type 2, the user must define EPS (2) and 
SIGMA (2) . 

For dispersive dielectric and magnetic materials, different 
material parameter arrays are used. The functional form of the 
frequency dependent permittivity/permeability that was 
implemented in the code is the Debye relaxation [4] with an 
effective DC conductivity given by 

e(co) = e'-je "= e c e + € 0 x e (<*>) + (16) 

3 <*> 


where the frequency dependent electric susceptibility function is 
defined as 
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x«(«> 


e s - e « 

1 + JWT 0 


(17) 


where e is the infinite frequency permittivity, e s is the zero 

<D 

frequency permittivity, r Q is the relaxation time, a is the 

effective electric conductivity, and to is the radian frequency. 
The same expressions were used for the frequency dependent 
permeability and can be written from equations (16) and (17) by 
replacing e by ix and o by o . The corresponding time domain 
susceptibility function is given by 


X e (t) 



(18) 


The FDTD implementation of frequency dependent permittivity 
and/or permeability involves a convolution with the electric 
and/or magnetic field and interested readers are referred to 
references [5-6] for further details. 

For dispersive dielectric materials, the corresponding 
material parameter arrays are EPSINF (e a ), EPSSTA (e s ) , RELAXT 

(t q ), and RELSIG (a). For dispersive magnetic materials, the 

arrays are XMUINF (Mj, XMUSTA (/i s ) , RELAXT ( t q ) and RELSIG (a). 

These dispersive material parameters are defined under the 
DISPERSIVE SETUP portion of the subroutine. The remainder of the 
subroutine computes constants used in field update equations and 
boundary conditions and writes the diagnostics file. 

SUBROUTINE EXSFLD 

This subroutine updates all x components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDONE array are used to determine the type of material present 
and the corresponding update equation to be used. These 
scattered field equations are based on the development given in 
[7]. 

SUBROUTINE EYSFLD 

This subroutine updates all y components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDTWO array are used to determine the type of material present 
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and the corresponding update equation to be used. 

SUBROUTINE EZSFLD 

This subroutine updates all z components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDTHRE array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINES RADEYX, RADEZX, RADEZY, RADEXY , RADEXZ and RADEYZ 

These subroutines apply the outer radiation boundary 
conditions to the scattered electric field on the outer 
boundaries of the problem space for non— magnetic scatterers. The 
user controls selection of these routines through the parameter 
MAGNET (described later) . 

SUBROUTINE HXSFLD 

This subroutine updates all x components of scattered 
magnetic field at each time step. IF statements based upon the 
IDFOR array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINE HYSFLD 

This subroutine updates all y components of scattered 
magnetic field at each time step. IF statements based upon the 
IDFIV array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINE HZSFLD 

This subroutine updates all z components of scattered 
magnetic field at each time step. IF statements based upon the 
IDSIX array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINES RADHXZ , RADHYX , RADHZY , RADHXY , RADHYZ and RADHZX 

These subroutines apply the outer radiation boundary 
conditions to the scattered magnetic field on the outer 
boundaries of the problem space for magnetic scatterers. The 
user controls selection of these routines through the parameter 
MAGNET (described later) . 


SUBROUTINE DATSAV * 
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This subroutine samples near zone scattered field quantities 
and saves them to disk. This subroutine is where the quantities 
to be sampled and their spatial locations are to be specified and 
is only called if near zone fields only are desired or if both 
near and far zone fields are desired. Total field quantities can 
also be sampled. See comments within the subroutine for 
specifying sampled scattered and/or total field quantities. When 
sampling magnetic fields, remember the S t/2 time difference 
between E and H when writing the fields to disk. Sections of 
code within this subroutine determine if the sampled quantities 
and the spatial locations have been properly defined. 

FUNCTIONS EXI, EYI , EZI, HXI, HYI and HZI 

These functions are called to compute the x, y and z 
components of incident electric and magnetic field, respectively. 
The functional form of the incident field is contained in a 
separate function SOURCE. 

FUNCTION SOURCE * 

This function contains the functional form of the incident 
field. The code as furnished uses the smooth cosine form of the 
incident field. An incident Gaussian pulse is also available by 
uncommenting the required lines and commenting out the smooth 
cosine pulse. Thus, this function need only be modified if the 
user changes the incident pulse from smooth cosine to Gaussian. 
Currently, only the smooth cosine pulse can be used for 
scattering from dispersive targets. A slight improvement in 
computing speed and vectorization may be achieved by moving this 
function inside each of the incident field functions EXI, EYI and 
so on. 

FUNCTIONS DEXI, DEYI , DEZI, DHXI , DHYI and DHZI 

These functions are called to compute the x, y and z 
components of the time derivative of incident electric and 
magnetic field, respectively. The functional form of the 
incident field is contained in a separate function DSRCE. 

FUNCTIONS DEXIXE , DEYIXE, DEZIXE, DHXIXE, DHYIXE and DHZIXE 

These functions compute the x, y and z components of the 
convolution of the time derivative of the incident field with the 

electric or magnetic Debye susceptibility function % • 


FUNCTION DSRCE * 



18 


This function contains the functional form of the time 
derivative of the incident field. The code as furnished uses the 
time derivative of the smooth cosine form of the incident field. 

A Gaussian pulse time derivative is also available by 
uncommenting the required lines and commenting out the smooth 
cosine pulse. Thus, the function need only be modified if the 
user changes from the smooth cosine to Gaussian pulse. Again, a 
slight improvement in computing speed and vectorization may be 
achieved by moving this function inside each of the time 
derivative incident field functions DEXI, DEYI and so on. 

FUNCTION DCONV 

This function evaluates the convolution of the time 
derivative of the incident field with the Debye susceptibility 

function % • 

SUBROUTINE ZERO 

This subroutine initializes various arrays and variables to 

zero. 

VIII. INCLUDE FILE DESCRIPTION (COMMOND. FOR) * 

The include file, COMMOND. FOR, contains all of the arrays 
and variables that are shared among the different subroutines. 
This file will require the most modifications when defining 
scattering problems. A description of the parameters that are 
normally modified follows. 

The parameters NX, NY and NZ specify the size of the problem 
space in cells in the x, y and z directions respectively. For 
problems where it is crucial to center the object within the 
problem space, then NX, NY and NZ should be odd for dielectric 
scatter ers and even for magnetic scatterers. The parameter NTEST 
defines the number of near zone quantities to be sampled and NZFZ 
defines the field output format. Set NZFZ=0 for near zone fields 
only, NZFZ=1 for far zone fields only and NZFZ=2 for both near 
and far zone fields. Parameter MAGNET is used to define magnetic 
scatterers and it controls the choice of RADE?? versus RADH?? 
absorbing boundary subroutines. It is set to 0 for dielectric 
scatterers and to 1 for magnetic scatterers. Parameter NUMMAT 
defines the total number of material types that are available for 
use. NEDISP and NHDISP define the number of dispersive 
dielectric and magnetic materials that are being used. Parameter 
NSTOP defines the maximum number of time steps. DELX, DELY, and 
DELZ (in meters) define the cell size in the x, y and z 
directions respectively. The 0 and <p incidence angles (in 
degrees) are defined by THINC and PHINC respectively and the 
polarization is defined by ETHINC and EPHINC. ETHINC=1.0, 
EPHINC=0.0 for 0-polarized incident field and ETHINC=0.0, 
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EPHINC=1.0 for 0-polarized incident fields. Parameters AMP and 
BETA define the maximum amplitude and the e temporal width of 
the incident pulse respectively. BETA automatically adjusts when 
the cell size is changed and normally should not be changed by 
the user. The far zone scattering angles are defined by THETFZ 
and PHIFZ . The code as furnished performs backscatter 
computations, but these parameters could be modified for a 
bistatic computation. 

IX. RCS COMPUTATIONS 

A companion code, RCS3D. FOR, has been included to compute 
RCS versus frequency. It uses the file name of the (FD) TD far 
zone output data (FZ0UT3D.DAT) and writes data files of far zone 
electric fields versus time (FZTIME.DAT) and RCS versus frequency 
(3DRCS.DAT) . The RCS computations are performed up to the 10 
cell/X 0 frequency limit. Refer to comments within this code for 
further details. 


X. RESULTS 


As previously mentioned, the code as furnished models a 6.67 
m radius, dispersive, Nickel Ferrite sphere and computes 
backscatter far zone scattered fields at angles of 0=22.5 and 
0=22.5 degrees. Results are included for 0.25 dB loaded 
dielectric and magnetic foam spheres, 60 dB dielectric and 
magnetic foam spheres and the Nickel Ferrite magnetic sphere. 

The material parameters for 0.25 dB and 60 dB loaded foam are: 


0.25 DB FOAM 


60 DB FOAM 


1.16 e s 

1.01 e B 


41.0 

1.6 


0.6497 ns r Q 

2 . 954E-04 S/m a 


0.3450 ns 
3 . 902E-01 S/m 


The magnetic materials are duals of the above as described 
earlier with the relative conductivity scaled by the ratio M 0 / e o- 

The Nickel Ferrite parameters are m b =1 / M s =100 r T o =20 ns and 

a*=0.0. For the 0.25 dB foams there are 10 cells per wavelength 
at approximately 3.0 GHz. For the 60 dB foam there are 10 cells 
per wavelength at approximately 2.35 GHz, as the 60 dB foam has a 
higher dielectric constant. 


Figures 2-10 show the real and imaginary parts of the 0.25 
dB and 60 dB foam permeability, the real and imaginary parts of 


the Nickel Ferrite permeability, and the magnetic 
susceptibilities versus time for all three materials. 
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Figures 11-14 show the co-polarized far zone electric field 
versus time and the co-polarized RCS for the 0.25 dB and 60 dB 
dielectric foam spheres respectively. 

Figures 15-18 show the co-polarized far zone electric field 
versus time and the co-polarized RCS for the 0.25 dB and 60 dB 
mangetic foam spheres respectively. 

Figures 19-20 show the co-polarized far zone electric field 
versus time and the co-polarized RCS for the default Nickel 
Ferrite sphere using the dispersive FDTD method and the non- 
dispersive FDTD method. For the non-dispersive method, an 
equivalent permeability and magnetic conductivity were defined at 
30 MHz from (16) with n replacing e as 

= ^ I (<a = 2*30E6) r a ^0 I (a>=2»30E6) 


XI. SAMPLE PROBLEM SETUP 

The code as furnished models a 6.67 m radius Nickel Ferrite 
dispersive magnetic sphere and computes backscatter far zone 
scattered fields at angles of 0=22.5 and 0= 22.5 degrees. The 
corresponding output data files are also provided, along with a 
code to compute Radar Cross Section using these data files. In 
order to change the code to a new problem, many different 
parameters need to be modified. A sample problem setup will now 
be discussed. 

Suppose that the problem to be studied is RCS backscatter 
versus frequency from a 28 cm by 31 cm perfectly conducting plate 
with a 3 cm dielectric coating with a dielectric constant of 4e 0 
using a 0-polarized field. The backscatter angles are 0=30.0 and 
0=60.0 degrees and the frequency range is up to 3 Ghz. 

Since the frequency range is up to 3 Ghz, the cell size must 
be chosen appropriately to resolve the field IN ANY MATERIAL at 
the highest frequency of interest. A general rule is that the 
cell size should be 1/10 of the wavelength at the highest 
frequency of interest. For difficult geometries, 1/20 of a 
wavelength may be necessary. The free space wavelength at 3 GHz 
is A 0 =10 cm and the wavelength in the dielectric coating at 3 GHz 
is 5 cm. The cell size is chosen as 1 cm, which provides a 
resolution of 5 cells/A in the dielectric coating and 10 cells/X 0 
in free space. Numerical studies have shown that choosing the 

cell size < 1/4 of the shortest wavelength in any material is 
the practical lower limit. Thus the cell size of 1 cm is barely 
adequate. The cell size in the x, y and z directions is set in 
the common file through variables DELX, DELY and DELZ. Next the 
problem space size must be large enough to accomodate the 
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scattering object, plus at least a five cell boundary (10 cells 
is more appropriate) on every side of the object to allow for the 
far zone field integration surface. It is advisable for plate 
scattering to have the plate centered in the x and y directions 
of the problem space in order to reduce the cross-polarized 
backscatter and to position the plate low in the z direction to 
allow strong specular reflections multiple encounters with the 
ORBC. A 10 cell border is chosen, and the problem space size is 
chosen as 49 by 52 by 49 cells in the x, y and z directions 
respectively. As an initial estimate, allow 2048 time steps so 
that energy trapped within the dielectric layer will radiate. 

Thus parameters NX, NY and NZ in COMMOND . FOR would be changed to 
reflect the new problem space size, and parameter NSTOP is 
changed to 2048. If all transients have not been dissipated 
after 2048 time steps, then NSTOP will have to be increased. 
Truncating the time record before all transients have dissipated 
will corrupt frequency domain results. Parameter NZFZ must be 
equal to 1 since we are interested in far zone fields only. 
Parameter MAGNET must be equal to 0 for the dielectric scatterer. 
To build the object, the following lines are inserted into the 
BUILD subroutine: 

C 

C BUILD THE DIELECTRIC SLAB FIRST 

C 

ISTART=1 1 

JSTART=11 

KSTART=11 

NXWIDE=28 

NYWIDE=3 1 

NZWIDE=3 

MTYPE=2 

CALL DCUBE ( ISTART , JSTART , K START , NXWIDE , NYWIDE , NZWIDE , MTYPE) 
C 

C BUILD PEC PLATE NEXT 

C 

I START= 1 1 

JSTART=11 

KSTART=11 

NXWIDE=28 

NYWIDE=3 1 

NZWIDE=0 

MTYPE=1 

CALL DCUBE (ISTART , JSTART , K START , NXWIDE , NYWIDE , NZWIDE , MTYPE ) 

The PEC plate is built last on the bottom of the dielectric 
slab to avoid any air gaps between the dielectric material and 
the PEC plate. In the common file, the incidence angles THINC 
and PHINC have to be changed to 30.0 and 60.0 respectively, the 
cell sizes (DELX , DELY, DELZ) are set to 0.01, and the 
polarization is set to ETHINC=1.0 and EPHINC=0.0 for 0-polarized 
fields. Since dielectric material 2 is being used for the 
dielectric coating, the constitutive parameters EPS (2) and 
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SIGMA ( 2 ) are set to 4e 0 and 0.0 respectively, in subroutine 
SETUP. This completes the code modifications for the sample 
problem. 

XII. NEW PROBLEM CHECKLIST 

This checklist provides a quick reference to determine if 
all parameters have been defined properly for a given scattering 
problem. A reminder when defining quantities within the code: 
use MKS units and specify all angles in degrees. 

COMMOND . FOR : 

1) Is the problem space sized correctly? (NX, NY, NZ) 

2) For near zone fields, is the number of sample points correct? 
(NTEST) 

3) Is parameter NZFZ defined correctly for desired field 
outputs? 

4) Is parameter MAGNET defined correctly for the type of 
scatterer? 

5) Is the number of dispersive dielectric (NEDISP) and 
dispersive magnetic (NHDISP) materials defined correctly? 

6) Is the number of time steps correct? (NSTOP) 

7) Are the cell dimensions ( DELX , DELY, DELZ) defined correctly? 

8) Are the incidence angles (THINC, PHINC) defined correctly? 

9) Is the polarization of the incident wave defined correctly 
(ETHINC, EPHINC)? 

10) For other than backscatter far zone field computations, are 
the scattering angles set correctly? (THETFZ, PHIFZ) 

SUBROUTINE BUILD: 

1) Is the object completely and correctly specified? 

SUBROUTINE SETUP: 

1) Are the constitutive parameters for each material specified 
correctly? (EPS, XMU, SIGMA, SIGMAC) 


2) Are the constitutive parameters for each dispersive material 
defined correctly? (EPSSTA, EPSINF, RELAXT , RELSIG, XMUINF, 
XMUSTA, RELAXT, RELSIG) 
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FUNCTIONS SOURCE and DSRCE: 

1) If the smooth cosine pulse is not desired, is it commented 

out and the Gaussian pulse uncommented? 

SUBROUTINE DATSAV : 

1) For near zone fields, are the sampled field types and spatial 

locations correct for each sampling point? (NTYPE, IOBS, JOBS, 

KOBS) 
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XI. FIGURE TITLES 

Fig. 1 Standard three dimensional Yee cell showing placement 
of electric and magnetic fields. 

Real part of relative permeability versus frequency for 
0.25 dB magnetic foam. 


Fig. 2 
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Fig. 3 Imaginary part of relative permeability versus 
frequency for 0.25 dB magnetic foam. 

Fig. 4 Relative magnetic susceptibility versus time for 0.25 
dB magnetic foam. 

Fig. 5 Real part of relative permeability versus frequency for 
60 dB magnetic foam. 

Fig. 6 Imaginary part of relative permeability versus 
frequency for 60 dB magnetic foam. 

Fig. 7 Relative magnetic susceptibility versus time for 60 dB 
magnetic foam. 

Fig. 8 Real part of relative permeability versus frequency for 
Nickel Ferrite with ju s = 100/ = 1/ r o = 20 ns * 

Fig. 9 Imaginary part of relative permeability versus 

frequency for Nickel Ferrite with /i s = 100, n m = 1, t q 

= 20 ns. 

Fig. 10 Relative magnetic susceptibility versus time for Nickel 
Ferrite. 

Fig. 11 Co-polarized far zone scattered field versus time for 
0.25 dB dielectric foam sphere with 20 cm radius. 

Fig. 12 Co-polarized Radar Cross Section versus frequency for 
0.25 dB dielectric foam sphere with 20 cm radius. 

Fig. 13 Co-polarized far zone scattered field versus time for 
0.25 dB magnetic foam sphere with 20 cm radius. 

Fig. 14 Co-polarized Radar Cross Section versus frequency for 
0.25 dB magnetic foam sphere with 20 cm radius. 

Fig. 15 Co-polarized far zone scattered field versus time for 
60 dB dielectric foam sphere with 20 cm radius. 

Fig. 16 Co-polarized Radar Cross Section versus frequency for 
60 dB dielectric foam sphere with 20 cm radius. 

Fig. 17 Co-polarized far zone scattered field versus time for 
60 dB magnetic foam sphere with 20 cm radius. 

Co-polarized Radar Cross Section versus frequency for 
60 dB magnetic foam sphere with 20 cm radius. 


Fig. 18 
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Fig. 19 Co-polarized far zone scattered field versus time for 
Nickel Ferrite sphere with 6.67 m radius using both 
dispersive FDTD and non-dispersive FDTD. 

Fig. 20 Co-polarized Radar Cross Section versus frequency for 
Nickel Ferrite sphere with 6.67 m radius using both 
dispersive FDTD and non-dispersive FDTD. 
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I . INTRODUCTION 

The Penn State Finite Difference Time Domain Electromagnetic 
Scattering Code Versions TEA and TMA are two dimensional 
numerical electromagnetic scattering codes based upon the Finite 
Difference Time Domain Technique (FDTD) first proposed by Yee [1] 
in 1966. The supplied version of the codes are two versions of 
our current two dimensional FDTD code set. This manual provides 
a description of the codes and corresponding results for the 
default scattering problem. The manual is organized into eleven 
sections: introduction, Version TEA and TMA code capabilities, a 

brief description of the default scattering geometry, a brief 
description of each subroutine, a description of the include 
files (TEACOM. FOR TMACOM. FOR) , a section briefly discussing 
scattering width computations, a section discussing the 
scattering results, a sample problem setup section, a new problem 
checklist, references and figure titles. 

II. VERSION TEA AND TMA CODE CAPABILITIES 

The Penn State University FDTD Electromagnetic Scattering 
Code Versions TEA and TMA have identical capabilities except the 
TEA code has electric field perpendicular to the z axis and the 
TMA code has electric field parallel to the z axis. Each code 
has the following capabilities: 

1) Ability to model lossy dielectric and perfectly conducting 
scatterers . 

2) First and second order outer radiation boundary condition 
(ORBC) operating on the electric fields for dielectric or 
perfectly conducting scatterers. 

3) Near to far zone transformation capability to obtain far zone 
scattered fields. 

4) Gaussian and smooth cosine incident waveforms with arbitrary 
incidence angles. 

5) Near zone field, current or power sampling capability. 

6) Companion codes for computing scattering width. 

III. DEFAULT SCATTERING GEOMETRY 

The codes as delivered are set up to calculate the far zone 
backscatter fields for an infinite, 0.25 m radius, perfectly 
conducting cylinder. The problem space size is 201 by 201 cells 
in the x and y directions, the cells are 1 cm squares, and the 
incident waveform is a Gaussian pulse with incidence angle of 
0=180 degrees. The output data files are included as a reference 
along with codes (SWTEA.FOR, SWTMA.FOR) for computing the 
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frequency domain scattering width using these output data files. 
The ORBC is the second order absorbing boundary condition set 
forth by Mur [2]. 

IV. SUBROUTINE DESCRIPTION 

In the description for each subroutine, an asterisk (*) will 
be placed by the subroutine name if that particular subroutine is 
normally modified when defining a scattering problem. Also, each 
subroutine will be denoted if it is applicable to the TE code 
only, the TM code only, or to both codes. 

MAIN ROUTINE (TE, TM) 

The main routine in the program contains the calls for all 
necessary subroutines to initialize the problem space and 
scattering object (s) and for the incident waveform, far zone 
transformation, field update subroutines, outer radiation 
boundary conditions and field sampling. 

The main routine begins with the include statement and then 
appropriate data files are opened, and subroutines ZERO, BUILD . 
and SETUP are called to initialize variables and/or arrays, build 
the object (s) and initialize the incident waveform and 
miscellaneous parameters, respectively. Subroutine SETFZ is 
called to intialize parameters for the near to far zone 
transformation if far zone fields are desired. 

The main loop is entered next, where all of the primary 
field computations and data saving takes place. During each time 
step cycle, the EXSFLD (TE) , EYSFLD (TE) , and EZSFLD (TM) 
subroutines are called to update the x, y, and z components of 
the scattered electric field. These scattered field equations 
are based upon the development given in [3]. RADEXY, RADEYX (TE) 
and RADEZX, RADEZY (TM) outer radiation boundary conditions are 
called next to absorb any outgoing scattered fields. Time is 
then advanced 1/2 time step according to the Yee algorithm and 
then the HXSFLD (TM) , HYSFLD (TM) , AND HZSFLD (TE) subroutines 
are called to update the x, y, and z components of scattered 
magnetic field. Time is then advanced another 1/2 step and then 
either near zone fields are sampled and written to disk in 
DATSAV , and/or the near zone to far zone vector potentials are 
updated in SAVFZ. The parameter NZFZ (described later) in the 
common file defines the type of output fields desired. 

After execution of all time steps in the main field update 
loop, subroutine FAROUT is called if far zone fields are desired 
to compute the far zone fields and write them to disk. At this 
point, the execution is complete. 
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SUBROUTINE SETFZ (TE, TM) 

This subroutine initializes the necessary parameters 
required for far zone field computations. The codes as furnished 
compute backscatter far zone field and can compute bistatic far 
zone fields for one scattering angle (i.e. one <p angle) . Refer 
to reference [4] for a complete description of the two 
dimensional near to far zone transformation. Other versions of 
this subroutine provide for multiple bistatic angles. 

SUBROUTINE SAVFZ (TE, TM) 

This subroutine updates the near zone to far zone vector 
potentials . 

SUBROUTINE FAROUT (TE, TM) 

This subroutine changes the near zone to far zone vector 
potentials to far zone electric field 0 (TM) and <p (TE) 
components and writes them to disk. 

SUBROUTINE BUILD (TE, TM) * 

This subroutine "builds" the scattering object (s) by 
initializing the IDONE (TE) , IDTWO (TE) , and IDTHRE (TM) arrays. 
The IDONE-IDTHRE arrays are for specifying perfectly conducting 
and lossy dielectric materials. Refer to Figure 1 for a diagram 
of the basic two dimensional Yee cell for the TE and TM case. 

For example, setting an element of the IDONE array at some I,J 
location is actually locating dielectric material at a cell edge 
whose center location is I+0.5,J. Thus, materials with diagonal 
permittivity tensors can be modeled. The default material type 
for all ID??? arrays is 0, or free space. By initializing these 
arrays to values other than 0, the user is defining an object by 
determining what material types are present at each spatial 
location. Other material types available for IDONE-IDTHRE are 1 
for perfectly conducting objects and 2-9 for lossy non-magnetic 
dielectrics. It is assumed throughout the code that all 
dielectric materials are non-magnetic (i.e. the materials have a 
permeability of /i 0 ) . This subroutine also has a section that 
checks the ID??? arrays to determine if legal material types have 
been defined throughout the problem space. The actual material 
parameters (e and a) are defined in subroutine SETUP. The 
default geometry is a 0.25 m radius, perfectly conducting 
cylinder . 

The user must be careful that his/her object created in the 
BUILD subroutine is properly formed. 

When it is important to place the object in the center of 
the problem space, NX and NY should be odd. This is due to the 



field locations in the Yee cell and also the placement of the E 
field absorbing boundary condition surfaces. 

If the object being modeled has curved surfaces, edges, etc. 
that are at an angle to one or more of the coordinate axes, then 
that shape must be approximately modeled by lines and faces in a 
"stair-stepped" (or stair-cased) fashion. This stair-cased 
approximation introduces errors into computations at higher 
freguencies. Intuitively, the error becomes smaller as more 
cells are used to stair-case a particular object. 

SUBROUTINE DPLATE (TE) 

This subroutine builds squares of dielectric material by 
defining two each of IDONE and IDTWO components corresponding to 
one spatial square of dielectric material. It can also be used 
to define thin (i.e. up to one cell thick) dielectric or 
perfectly conducting wires. Refer to comments within DPLATE for 
a description of the arguments and usage of the subroutine. 

SUBROUTINE SETUP (TE, TM) * 

This subroutine initializes many of the constants required 
for incident field definition, field update equations, outer 
radiation boundary conditions and material parameters. The 
material parameters e and a are defined for each material type 
using the material arrays EPS and SIGMA respectively. The array 
EPS is used for the total permittivity and SIGMA is used for the 
electric conductivity. These arrays are initialized in SETUP to 
free space material parameters for all material types and then 
the user is required to modify these arrays for his/her 
scattering materials. Thus, for the lossy dielectric material 
type 2, the user must define EPS(2) and SIGMA (2) . The remainder 
of the subroutine computes constants used in field update 
equations and boundary conditions and writes the diagnostics 
file. 

SUBROUTINE EXSFLD (TE) 

This subroutine updates all x components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDONE array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINE EYSFLD (TE) 

This subroutine updates all y components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDTWO array are used to determine the type of material present 
and the corresponding update equation to be used. 
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SUBROUTINE EZSFLD (TM) 

This subroutine updates all z components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDTHRE array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINES RADEYX, RADEXY (TE) and RADEZX, RADEZY (TM) 

These subroutines apply the outer radiation boundary 
conditions to the scattered electric field on the outer 
boundaries of the problem space. 

SUBROUTINE HXSFLD (TM) 

This subroutine updates all x components of scattered 
magnetic field at each time step. The standard non-magnetic 
update equation is used. 

SUBROUTINE HYSFLD (TM) 

This subroutine updates all y components of scattered 
magnetic field at each time step. The standard non-magnetic 
update equation is used. 

SUBROUTINE HZSFLD (TE) 

This subroutine updates all z components of scattered 
magnetic field at each time step. The standard non-magnetic 
update equation is used. 

SUBROUTINE DATSAV (TE, TM) * 

This subroutine samples near zone scattered field quantities 
and saves them to disk. This subroutine is where the quantities 
to be sampled and their spatial locations are to be specified and 
is only called if near zone fields only are desired or if both 
near and far zone fields are desired. Total field quantities can 
also be sampled. See comments within the subroutine for 
specifying sampled scattered and/or total field quantities. When 
sampling magnetic fields, remember the <St/2 time difference 
between E and H when writing the fields to disk. Sections of 
code within this subroutine determine if the sampled quantities 
and the spatial locations have been properly defined. 

FUNCTIONS EXI, EYI (TE) and EZI (TM) 


These functions are called to compute the x, y and z 
components of incident electric field. The functional form of 
the incident field is contained in a separate function SOURCE. 



8 


FUNCTION SOURCE (TE, TM) * 

This function contains the functional form of the incident 
field. The code as furnished uses the Gaussian form of the 
incident field. An incident smooth cosine pulse is also 
available by uncommenting the required lines and commenting out 
the Gaussian pulse. Thus, this function need only be modified if 
the user changes the incident pulse from Gaussian to smooth 
cosine. A slight improvement in computing speed and 
vectorization may be achieved by moving this function inside each 
of the incident field functions EXI, EYI and so on. 

FUNCTIONS DEXI, DEYI (TE) and DEZI (TM) 

These functions are called to compute the x, y and z 
components of the time derivative of incident electric field. 

The functional form of the incident field is contained in a 
separate function DSRCE. 

FUNCTION DSRCE (TE, TM) * 

This function contains the functional form of the time 
derivative of the incident field. The code as furnished uses the 
time derivative of the Gaussian form of the incident field. A 
smooth cosine pulse time derivative is also available by 
uncommenting the required lines and commenting out the Gaussian 
pulse. Thus, the function need only be modified if the user 
changes from the Gaussian to smooth cosine pulse. Again, a 
slight improvement in computing speed and vectorization may be 
achieved by moving this function inside each of the time 
derivative incident field functions DEXI, DEYI and so on. 

SUBROUTINE ZERO (TE, TM) 

This subroutine initializes various arrays and variables to 

zero. 

V. INCLUDE FILE DESCRIPTION ( TEACOM . FOR , TMACOM . FOR ) * 

The include files, TEACOM. FOR, TMACOM. FOR, contain all of 
the arrays and variables that are shared among the different 
subroutines. These files will require the most modifications 
when defining scattering problems. A description of the 
parameters that are normally modified follows. 

The parameters NX and NY specify the size of the problem 
space in cells in the x and y directions respectively. For 
problems where it is crucial to center the object within the 
problem space, then NX and NY should be odd. The parameter NTEST 
defines the number of near zone quantities to be sampled and NZFZ 
defines the field output format. Set NZFZ=0 for near zone fields 
only, NZFZ=1 for far zone fields only and NZFZ=2 for both near 
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and far zone fields. Parameter NSTOP defines the maximum number 
of time steps. DELX and DELY (in meters) define the cell size in 
the x and y directions respectively. The 0 incidence angle ( in 
degrees) is defined by PHINC and the polarization is defined by 
the code that is being used. Parameters AMP and BETA define the 
maximum amplitude and the e temporal width of the incident 
pulse respectively. BETA automatically adjusts when the cell 
size is changed and normally should not be changed by the user. 
The far zone scattering angle is defined by PHIFZ. The codes as 
furnished perform backscatter computations, but this parameter 
could be modified for a bistatic computation. 

VI. SCATTERING WIDTH COMPUTATIONS 

Companion codes, SWTEA . FOR , SWTMA . FOR , have been included to 
compute scattering width versus freguency. Each code uses the 
file name of the FDTD far zone output data (FZOUTTE.DAT, 
FZOUTTM.DAT) and writes data files of far zone electric field 
versus time (FZTE.DAT, FZTM.DAT) and scattering width versus 
frequency (SWTE.DAT, SWTM.DAT). The scattering width 
computations are performed up to the 10 cell/A 0 frequency limit. 
Refer to comments within these codes for further details. 

VII . RESULTS 

As previously mentioned, the codes as furnished model an 
infinite, 0.25 m radius, perfectly conducting cylinder and 
compute backscatter far zone scattered field at an angles of 
0=180 degrees. 

Figures 2-4 show the TM far zone electric field versus time 
and the scattering width magnitude and phase for the 0.25 m 
radius perfectly conducting cylinder. Keep in mind that the far 
zone time domain electric field shown here is not the actual time 
domain scattered field. The actual far zone time domain 
scattered field can be obtained by an FFT of the FDTD time domain 
results, then multiplying the FFT by the appropriate frequency 
domain factor (described in [4]) and performing an inverse FFT. 

Figures 5-7 show the TE far zone electric field versus time 
and the scattering width magnitude and phase for the 0.25 m 
radius perfectly conducting cylinder. 

VIII. SAMPLE PROBLEM SETUP 

The codes as furnished model an infinite, 0.25 m radius, 
perfectly conducting cylinder and compute backscatter far zone 
scattered field at an angle of 0=180 degrees. The corresponding 
output data files are also provided, along with codes to compute 
scattering width using these data files. In order to change the 
code to a new problem, many different parameters need to be 
modified. A sample problem setup will now be discussed. 
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Suppose that the problem to be studied is scattering width 
backscatter versus frequency from a 0.30 m radius dielectric 
cylinder with a dielectric constant of 4e 0 using a 0-polarized 
field. The backscatter angle is 0=60.0 degrees and the frequency 
range is up to 3 Ghz . 

The incident field is 0-polarized which indicates the TM 
code must be used. Since the frequency range is up to 3 Ghz, the 
cell size must be chosen appropriately to resolve the field IN 
ANY MATERIAL at the highest frequency of interest. A general 
rule is that the cell size should be 1/10 of the wavelength at 
the highest frequency of interest. For difficult geometries, 

1/20 of a wavelength may be necessary. The free space wavelength 
at 3 GHz is A 0 =10 cm and the wavelength in the dielectric coating 
at 3 GHz is 5 cm. The cell size is chosen as 1 cm, which 
provides a resolution of 5 cells/ A. in the dielectric coating and 
10 cells/ A 0 in free space. Numerical studies have shown that 

choosing the cell size < 1/4 of the shortest wavelength in any 
material is the practical lower limit. Thus the cell size of 1 
cm is barely adequate. The cell size in the x and y directions 
is set in the common file through variables DELX and DELY. Next 
the problem space size must be large enough to accomodate the 
scattering object, plus at least a five cell boundary (10 cells 
is more appropriate) on every side of the object to allow for the 
far zone field integration surface. The default problem space 
size of 201 by 201 is adequate and provides a 75 cell border 
around the cylinder. As an initial estimate, allow 2048 time 
steps so that energy trapped within the dielectric layer will 
radiate. Thus parameter NSTOP is changed to 2048. If all 
transients have not been dissipated after 2048 time steps, then 
NSTOP will have to be increased. Truncating the time record 
before all transients have dissipated will corrupt frequency 
domain results. Parameter NZFZ must be equal to 1 since we are 
interested in far zone fields only. To build the object, simply 
change the RADIUS variable in the BUILD subroutine to 0.30. In 
the common file, the incidence angle PHINC has to be changed to 
60.0 respectively, and the cell sizes (DELX and DELY) are set to 
0.01. Since dielectric material 2 is being used for the 
dielectric coating, the constitutive parameters EPS (2) and 
SIGMA ( 2 ) are set to 4e 0 and 0.0 respectively, in subroutine 
SETUP. This completes the code modifications for the sample 
problem. 

IX. NEW PROBLEM CHECKLIST 

This checklist provides a quick reference to determine if 
all parameters have been defined properly for a given scattering 
problem. A reminder when defining quantities within the code: 
use MKS units and specify all angles in degrees. 



TEACOM . FOR , TMACOM . FOR : 

1) Is the problem space sized correctly? (NX, NY) 
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2) For near zone fields, is the number of sample points correct? 
(NTEST) 

3) Is parameter NZFZ defined correctly for desired field 
outputs? 

4) Is the number of time steps correct? (NSTOP) 

5) Are the cell dimensions ( DELX , DELY) defined correctly? 

6) Is the incidence angle (PHINC) defined correctly? 

7) For other than backscatter far zone field computations, is 
the scattering angle set correctly? (PHIFZ) 

SUBROUTINE BUILD: 

1) Is the object completely and correctly specified? 

SUBROUTINE SETUP: 

1) Are the constitutive parameters for each material specified 
correctly? (EPS and SIGMA) 

FUNCTIONS SOURCE and DSRCE : 

1) If the Gaussian pulse is not desired, is it commented out and 
the smooth cosine pulse uncommented? 

SUBROUTINE DATSAV : 

1) For near zone fields, are the sampled field types and spatial 
locations correct for each sampling point? (NTYPE, IOBS, JOBS) 
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XI. FIGURE TITLES 

Fig. 1 Standard two dimensional Yee cell showing placement of 
electric and magnetic fields for the TE and TM case. 

Fig. 2 Far zone scattered field versus time for 0.25 m radius 
perfectly conducting cylinder using TM polarization. 

Fig. 3 Scattering width magnitude versus frequency for 0.25 m 
radius perfectly conducting cylinder using TM 
polarization. 

Fig. 4 Scattering width phase versus frequency for 0.25 m 
radius perfectly conducting cylinder using TM 
polarization. 

Fig. 5 Far zone scattered field versus time for 0.25 m radius 
perfectly conducting cylinder using TE polarization. 

Fig. 6 Scattering width magnitude versus frequency for 0.25 m 
radius perfectly conducting cylinder using TE 
polarization. 

Fig. 7 Scattering width phase versus frequency for 0.25 m 
radius perfectly conducting cylinder using TE 
polarization. 
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Abstract 


n a previous paper [1] a time domain transformation useful 
or extrapolating three dimensional near zone finite difference 
time domain (FDTD) results to the far zone was presented. In 
this paper the corresponding two dimensional transform is 
outlined. While the three dimensional transformation produced a 
physically observable far zone time domain field, this is not 
convenient to do directly in two dimensions, since a convolution 
would be required. However, a representative two dimensional far 
zone time domain result can be obtained directly. This result 
can then be transformed to the frequency domain using a Fast 

r ^ er ^ rans ^ orin / corrected with a simple multiplicative factor 
and used, for example, to calculate the complex wideband 
scattering width of a target. if an actual time domain far 
result is required it can be obtained by inverse Fourier 
transform of the final frequency domain result. 


zone 



Introduce j nn 


A previous paper [ 1 ] described a method for transforming 
near zone Finite Difference Time Domain (FDTD) results directly 

° far ZOne without «tst transforming to the frequency 
omain. This far zone field could then be used to compute the 
scattering cross section of an illuminated target or an antenna 
radiation pattern over the entire frequency band of the FDTD 
calculations. A similar result was published in [2], m this 
paper the corresponding two dimensional transform is presented 
and validated by comparison with calculated results for a 
perfectly conducting circular cylinder. 


Approach 


in [1] the frequency domain far zone transformation 

IZTZl H FOUriSr tranSf ° raed to the doma in and used in 

O erive an approach to transform near zone FDTD 

fields to the far zone directly in the time domain. Our approach 

ere will be to present the fundamental frequency domain 

equations for both two and three dimensions, and by comparing 

them obtain the factor needed to convert the three dimensional 

far zone transform to function in two dimensions. 

We again surround the scatterer with a closed surface s' 

and consider that equivalent tangential electric and magnetic' 

time harmonic surface currents may exist on this surface 

Referring to f 3], we obtain the vector potentials for the three 
dimensional case as 
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-jkr 


( 2 ) 


4nr 


J M s e jkfi ' f ds J 

s' 


with 



k 


= 0 ) 



(the wave number) , f the unit vector 


to the far zone field point, r / the vector to the source point 

of integration, r the distance to the far zone field point, and 
S' the closed surface surrounding the scatterer. 

The far zone frequency domain electric fields of the 
scatterer are then obtained from 

E e = Ag - jw //ie" F 0 (3) 


E * = - jw n + 



(4) 


One can then easily convert to radar cross section (RCS) , if 
desired, by applying 


3D 
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(47rr 


S , 2 

E 3 dI 
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( 5 ) 


where Ejd is either E 0 or E 0 of (2) or (3), and E 1 ’ is the 
incident plane wave electric field. 


In [1] the equations corresponding to (1,2) were easily 
transformed to the time domain since the exponential phase term 
inside the integrals corresponds to a time shift relative to an 
arbitrary time reference point. Equations corresponding to (2,3) 
were also readily transformed to the time domain since the jo 
factor in these equations corresponds to a time derivative. Thus 
ulting ^ lne domain fields can be obtained conveniently 
directly in the time domain, as shown in [1], 


Now consider the corresponding equations 
The vector potentials are given by 


in two dimensions. 
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e -jkp 
^8 jk7rp 


J j $ e ikp «■<#-♦'> ds ' 
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here p and <p are the coordinates of the source point of 

integration, and p and 0 the coordinates of the far zone field 

point. The corresponding far zone radiated fields are obtained 
from 


E z - -jo m A z + jo> \j F 


(8) 


E * - - A - jo {pie F 


( 9 ) 


Also, the two dimensional scattering width is defined as 

lim 


a — ~ , l^l 2 

= p- (2 " P 7^7I 


) 


( 10 ) 



where Ejo 


is either E z or E t of (8,9). 


The approach applied in [l] cannot be conveniently applied 
in the two dimensional case due to the factor of 1/y/JTc 

(actually 1/VjWSe ) in (6,7). in order to evaluate the 

Fourier transform of (6,7) directly in the time domain a 
convolution operation would be required. To avoid this 
complication our approach will be to modify the results in [ 1 ] to 
provide representative two dimensional time domain far zone 
fields which can then be converted to the actual frequency domain 
fields by a multiplication in the frequency domain rather than 
the time domain convolution. Should the actual time domain far 
zone fields be required, they can then be obtained by an 

additional Fourier transform of these results back to the time 
domain. 

In order to convert our previous three dimensional results 
to two dimensions, we compare the two sets of equations. First, 
comparing (3,4) with (8,9), since the spherical unit 0 vector is 

equal to the negative of the cylindrical unit z vector, (3,4) and 

(8,9) correspond exactly, and no adjustment between two and three 
dimensional transforms is needed. 

Next, comparing (1,2) with (6,7) the r 2 and p factors are 
compensated by the definitions in (5) and (10) respectively, as 
expected, and no compensation is needed here either. 

Finally, consider equations (1,2) vs (6,7). The additional 
dimension of integration in (1,2) is compensated for by defining 
a scattering width per unit length (in z) in (10). This 


corresponds in (1,2) to having no z variation and integrating the 
2 variable over a unit distance. The exponents provide 
equivalent phase (time) delays and need not be compensated for. 
Considering the remaining factors, it is easily determined that 
in the frequency domain, the relationship between far zone 
electric fields obtained from a three dimensional far zone 
transformation with no z variation and the two dimensional far 
zone fields is 
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2nc 

jw 



( 11 ) 


where 


C = 1/VM6 


With these results the time domain far zone transform given 
in [l] can be easily adapted to two dimensions as follows: 

1) Consider only the field components and corresponding surface 
currents excited in the two dimensional problem. For example, 
for a TE Z computation only H z , E x , E y , and the corresponding 
surface currents are included. 

2) Calculate the far zone time domain fields using the method 
described in [1], but for a two dimensional integration surface 
which encloses the scatterer. Let Sz, the z coordinate unit cell 
dimension used in [1], equal l (meter). (This field is not a 
physically observable field. It represents the radiation from a 

l en 9 th of the scatterer in the time domain.) 

3) Fourier transform the result of step 2) and multiply the 
result by the factor in (11) . This result is the frequency 
domain two dimensional far zone field, which can then be used in 



(10) to calculate the scattering width as a function of 
frequency . 

4) if the actual time domain two dimensional far zone field is 
desired, it can be obtained by an additional Fourier 
transformation of the result obtained in 3) back to the time 
domain. 

Demonstration 

In order to demonstrate the capabilities of the above 
approach a pair (TE Z and TM Z ) of FDTD codes were developed from 
the three dimensional FDTD code described in [l]. These codes 
utilize second order Mur absorbing boundary conditions acting on 
the electric fields. The test geometry was a circular perfectly 
conducting cylinder of radius 0.25 meters. 

Both TE and TM polarization was considered, and two sets of 
calculations were made for each polarization. For the first set 
the FDTD cells were 1 cm squares, with a problem space 200x200 
cells. For the second set the cells were 0.5 cm squares in a 
500x500 cell problem space. On a 25 MHz 486 PC (approximately 1 
MFLOP) each of the first set required about 20 minutes to 
compute, with each of the second set requiring a few minutes less 
than 2 hours . 

For all cases the incident plane wave traveled in the x 
direction, backscatter was calculated, and 2048 time steps were 
evaluated. In order to clearly show the response, not all time 
steps are included in the Figures. 

Figure 1 shows the relative far zone fields computed 
directly in the time domain as outlined above for the TM 
polarization. The small ripple at approximately 15 ns on Figure 
1 is due to reflections from the Mur outer boundary. Figures 2 
and 3 show comparison with the exact solution for the scattering 


width amplitude and phase. The upper frequency limit of 3.0 GHz 
corresponds to 10 cells per wavelength. The agreement is quite 
good. 


A similar set of data for TE polarization is shown in 
Figures 4-6. The time domain far zone results in Figure 4 has 
ripples in the 6-8 ns range due to the staircasing of the round 
cylinder with square FDTD cells. The small negative pulse at 10 
ns is the creeping wave which has traveled around the cylinder. 
Again, there is a small ripple at 14 ns due to the Mur boundary 
reflection. This is the difficult polarization for approximating 
a smooth surface with a "staircased" FDTD code, yet the agreement 
in Figures 5 and 6 is reasonably good, reproducing the first 6 
1/2 ripples in the scattering width. 

The above calculations are repeated in Figures 7-13 with a 
greater expenditure of computer resources. For these results the 
cell size of 0.5 cm changes the 3.0 GHz upper frequency limit of 
the plots to correspond to 20 cells per wavelength. The 
improvement in the agreement with the exact solution is clear, 
indicating the accuracy that can be obtained from this approach. 
Note especially Figure 12, which shows on a expanded dB scale 
agreement with the exact solution within a fraction of a dB for 
the first 9 lobes of the response. 


Conclusions 

A simple approach to calculating a wide bandwidth time 
domain transformation of near zone FDTD fields to the far zone 
has been presented. It is based on simple modifications to the 
previous three dimensional method presented in [1]. Results 
obtained using this transformation show good agreement with the 
exact solution for a circular cylinder for both polarizations. 
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Far Zone relative electric field vs time for TM Z 
polarized incident Gaussian pulsed plane wave 
illuminating a perfectly conducting circular cylinder. 
FDTD cells are l cm squares. 

Scattering width amplitude obtained from far zone time 
domain results and compared with exact solution. FDTD 
cells are l cm squares. 

Scattering width phase obtained from far zone time 
domain results and compared with exact solution. FDTD 
cells are l cm squares. 

Far Zone relative electric field vs time for TE Z 
polarized incident Gaussian pulsed plane wave 
illuminating a perfectly conducting circular cylinder. 
FDTD cells are 1 cm squares. 

Scattering width amplitude obtained from far zone time 
domain results and compared with exact solution. FDTD 
cells are 1 cm squares. 

Scattering width phase obtained from far zone time 
domain results and compared with exact solution. FDTD 
cells are l cm squares. 

Far Zone relative electric field vs time for TM Z 
polarized incident Gaussian pulsed plane wave 
illuminating a perfectly conducting circular cylinder. 
FDTD cells are 0.5 cm squares. 

Scattering width amplitude obtained from far zone time 
domain results and compared with exact solution. FDTD 
cells are 0.5 cm squares. 



Fig 9. Scattering width phase obtained from far zone time 

domain results and compared with exact solution. FDTD 
cells are 0.5 cm squares. 

Fig 10 . Far Zone relative electric field vs time for TE Z 
polarized incident Gaussian pulsed plane wave 
illuminating a perfectly conducting circular cylinder. 
FDTD cells are 0.5 cm squares. 

Fig li. Scattering width amplitude obtained from far zone time 
domain results and compared with exact solution. FDTD 
cells are 0.5 cm squares. 

Fig 12. Scattering width amplitude obtained from far zone time 
domain results and compared with exact solution on an 
expanded dB scale. FDTD cells are 0.5 cm squares. 

13. Scattering width phase obtained from far zone time 

domain results and compared with exact solution. FDTD 
cells are 0.5 cm squares. 
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ABSTRACT 

Surface impedance boundary conditions are employed to reduce 
the solution volume during the analysis of scattering from lossy 
dielectric objects. In a finite difference solution, they also can 
be utilized to avoid using small cells, made necessary by shorter 
wavelengths in conducting media throughout the solution volume. 
The standard approach is to approximate the surface impedance over 
a very small bandwidth by its value at the center frequency, and 
then use that result in the boundary condition. In this paper, two 
implementations of the surface impedance boundary condition are 
presented. One implementation is a constant surface impedance 
boundary condition and the other is a dispersive surface impedance 
boundary condition that is applicable over a very large frequency 
bandwidth and over a large range of conductivities. Frequency 
domain results are presented in one dimension for two conductivity 
values and are compared with exact results. Scattering width 
results from an infinite square cylinder are presented as a two 
dimensional demonstration. Extensions to three dimensions should 
be straightforward. 

I. Introduction 

The Finite Difference Time Domain (FDTD) technique permits the 
analysis of interactions of electromagnetic waves with objects of 
arbitrary shape and material composition. This method was first 
proposed by Yee [1] for isotropic, non-dispersive materials in 
1966; and through various modifications during the past twenty 
years, it has evolved into a mature computational technique. 
Reference [2] and the references contained therein provide an 
account of various extensions and modifications of the original 
FDTD algorithm. The present FDTD technique is capable of 
electromagnetic scattering analysis from objects of arbitrary and 
complicated geometrical shape and material composition over a large 
band of frequencies. This technique has recently been extended to 
include dispersive dielectric materials [3], chiral materials [4] 
and plasmas [5]. Due to these numerous capabilities, the FDTD 
method has begun to gain widespread acceptance as a viable 
computational alternative to the classical method of moments (MM) 
technique for many problems. 

To analyze electromagnetic field interaction with lossy 
dielectric objects, the FDTD method requires that the interior of 
the object be modeled in order for fields to penetrate the body. 
Accurate modeling often requires a very fine spatial grid resulting 
in a relatively large number of cells for moderately sized objects. 
A highly conducting dielectric object can be replaced by a surface 
impedance boundary condition (SIBC) that is a function of the 
material parameters. Thus, this boundary condition eliminates the 
spatial quantization of the object and reduces the overall size of 
the solution space not only by eliminating cells within the lossy 
dielectric, but also by allowing larger cells to be used in the 
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exterior region* As with any computational electromagnetic tool, a 
technique that reduces the solution space or number of unknowns is 
quite welcome. 

Of historical interest, surface impedance boundary conditions 
were first proposed by Leontovich in the 1940’s [6] and were 
rigorously developed by Senior in 1960 [7]. During the past thirty 
years, researchers have applied surface impedance concepts in the 
frequency domain to numerous electromagnetic scattering problems. 
Time domain surface impedance concepts received little attention 
until recently. Through some impressive work, Maloney and Smith 
[8] have previously implemented a surface impedance boundary 
condition in the FDTD method. However, their implementation has a 
minor disadvantage because the exponential rates and coefficients 
for recursive updating have to be reevaluated each time the 
conductivity or loss tangent is changed. With our proposed method, 
the exponential rates and coefficients only have to be evaluated 
once . Tesche [9] has also investigated surface impedance concepts 
in an integral equation time domain solution, but presented limited 
computational time domain results. 

It is the purpose of this paper to introduce a constant 
surface impedance boundary condition that is applicable for a 
single frequency and a dispersive surface impedance boundary 
condition that is applicable over a large frequency bandwidth and 
range of conductivities. The dispersive surface impedance includes 
frequency variations which results in a time domain boundary 
condition involving a convolution. We will then show how to 
efficiently evaluate this convolved surface impedance using 
recursion. 

II. Motivation 

The motivation for implementing a SIBC in the FDTD method is 
to reduce the computational resource requirements for modeling 
highly conducting lossy dielectric objects. In the standard FDTD 
method, modeling highly conducting lossy dielectric objects 
requires that the cell size be chosen small enough to resolve the 
field inside the object at the maximum frequency of interest. For 
example, suppose scattering from a lossy dielectric object with 
permeability /i, permittivity € and conductivity a=2 . 0 S/m is to be 
studied over the frequency band 0-10 GHz. The cell size must be 
chosen as some fraction of the wavelength inside the conducting 
material at the maximum frequency of interest. Thus the cell size 
is chosen (typically) as 


«5x = 6y = 5z = — 
10 


( 1 ) 


4 


where £ r is the complex relative permittivity constant of the 

material and A. and A 0 are the wavelengths inside the material and 
in free space at 10 GHz, respectively. The complex permittivity 
for lossy dielectrics in the frequency domain is 

a = e + — (2) 

ju 


where w is the radian frequency. The complex relative permittivity 
is determined using (2) as 



£ r + 


a 

j^o 


(3) . 


If the material is a good conductor over all frequencies of 
interest, then the constitutive parameters satisfy the condition 

— » 1 (4) . 

(i)€ 


Therefore, 2 r can be approximated as 


Assuming parameters M = Mo anc * 6=e o anc * us i n 9 the values of £ r and A 0 

at 10 GHz, the cell size is Sx = Sy = Sz = 1.582 mm. If a SIBC is 
used, then the cell size need only be chosen to resolve the field 
in free space and (1) is modified to 

^ x sibc = 5 Ysibc = ^ z sibc _ 


Again, using the value for 1 0 at 10 GHz, the cell size is £x S]BC - 
$y slB c = <Sz slBC = 3.0 mm. Thus the cell size has been increased by 

the factor ^ | =1.90, and the computational storage requirements 

are reduced by the same factor. Therefore, the computational 
savings, denoted by S, is 


S 



(V) 


where e r is given by (5), (4) is satisfied for all frequencies of 

interest and d is the number of dimensions. 
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III. FDTD Constant Surface Impedance Implementation 

To implement the constant SIBC in the FDTD method we consider 
the planar air-lossy dielectric interface as shown in Figure 1. 
The conducting material has permittivity e, permeability y, and 
conductivity cr. We assume that the thickness of the material is 
large compared to the skin depth. We will also assume that the 
material is linear and isotropic and a basic familiarity with the 
Yee algorithm [1]. Figure 1 also shows the one-dimensional FDTD 
grid. 

The first order (or Leontovich) impedance boundary condition 
relates tangential total field components and is given in the 
frequency domain as [6] 

E x (o>) = Z s (a>)H y (o>) (8) 


where Z (<a) is the surface impedance of the conductor, 
frequency domain surface impedance for good conductors is 


Z s (o) 


(l + j) 


o)j a 


yi 2a 




N 


(9) • 


Using (9) , (8) can be rewritten as 

E x ( g>) =(R s (<i>) + j X s ( o>) ) H y ( o>) 


( 10 ) 


where R s is the surface resistance and X s is the surface reactance. 
Consider rewriting (10) as 


E x (o» 



j<i>L s (o>) j 


H y (o>) 


(ID 


with the resistance and inductance defined by 


R s (g)) 


(O/i 


nJ 2a 


M 


( 12 ) 


L s (g>) 


'sj 2 aoi 
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To remove the frequency dependence of the surface resistance and 
inductance, these quantities are evaluated at a particular 
frequency and are subsequently treated as constants. Equation (11) 
then becomes 

E x (q) = (r s + jo>L s )H y (o) (13) . 


This is the required frequency domain constant surface impedance 
boundary condition. To incorporate this boundary condition into 
the FDTD algorithm, the time domain equivalent of (13) must be 
obtained. Performing an inverse Fourier transform operation on 
(13) results in 


E,(t) 


R,H y (t) +L s^ H v (t) 


(14) . 


This equation defines the time domain FDTD constant surface 
impedance boundary condition. 

To implement this constant surface impedance boundary 
condition, space and time are quantized by defining 

z -* (k<Sz) =» (k) (15) # 

t - (n<St) =* (n) 


The Faraday-Maxwell law is then used to obtain the H y component in 
the free space cell next to the impedance boundary. Since the 
impedance boundary condition requires that the electric and 
magnetic fields are co-located in space and time, we assume that 
the magnetic field 1/2 cell in front of the impedance boundary and 
1/2 time step previous is an adequate approximation. The Faraday- 
Maxwell law yields 


-(m 0 $x<S z) 


d(nSt) 


HY n (k + 1/2 ) 


EX n (k + l)5x - EX n (k) Sx (16) 


Note the component EX n (k+l) of (16) is the electric field component 
at the impedance boundary. Quantizing space and time in (14) and 
using the result to eliminate EX n (k+l) in (16) gives 


-(m 0 *z) 


a (nat) 


HY n (k + 1/ 2 ) 


= R HY n (k + 1/2 ) + L 


8(n£t) 


HY n (k + 1/2 ) 


- EX n (k) 
(17) • 


Notice that the HY n (k+1/2) term in (17) is time indexed at time 
step n. This term is approximated as 
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HY n (k + 1/ 2 ) « - HY n+1/2 (k + 1/2 ) + HY n ' 1/2 (k + 1/ 2 ) J (18). 

2 

Using (18) , and approximating the time derivatives on the magnetic 
fields in (17) as finite differences gives 

-(m 0 5z +L s )(HY n+1/2 (k + l/2) -HY n ' 1/2 (k + 1/2 )) = -y- (hy" + 1/2 (k + 1/ 2 ) - HY n ~ 1/2 (k + 1/2 ) ) 

- StEX n (k) (19) . 

Solving for HY n+V2 (k+l/2) in (19) yields 


HY n+1/2 (k + 1/2) 


H n Sz +L, -R ( .6t/2| 

= ! ! HY n ' /2 (k + 1/2 ) 

H 0 6z +L S +R s <St/2 


St 

M 0 5z+L s +R s 6t/2 


EX n (k) 
( 20 ) . 


This equation implements the constant surface impedance boundary 
condition in the FDTD method. 

IV. FDTD Dispersive Surface Impedance Implementation 

To derive a similar relation to (20) valid over a wide 
frequency band, we begin with the same set of underlying 
assumptions as for the constant surface impedance. The primary 
exception is that the surface impedance will vary with frequency 
and will not be approximated by its value at a particular 
frequency. All frequency domain information is inversed Fourier 
transformed to equivalent time domain form. The SIBC is then 
implemented in the FDTD method with the required convolution using 
a recursive updating technique. 

The standard first order impedance boundary condition remains 
unchanged and is given by (8). In a similar fashion as Tesche [9], 
(8) is rewritten as 

Z s (<o) 

E ( (j>) = jo> H (o) 

jco J 


( 21 ) . 
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Defining 


Zs(<*>) 


z s (“> 

jw 


and substituting (9) into (22) gives 



( 22 ) 


(23) . 


Substituting (23) into (21), a modified surface impedance boundary 
condition is obtained as 

E x (u>) = z'(w) [j«H y (<o) ] ( 24 ) • 


The time domain equivalent of (24) is obtained via an inverse 
Fourier transform operation as 

E x (t) = Zg ( t ) * -A H y (t) (25) 

[ dt 

where the asterisk denotes convolution, 

E x (t) = ^•- 1 [E x (o>) ] 

H y (t) = & [ H y ( <i>) ] (26) 

Zg(t) = ^^[Zgto)) ] 


and the denotes the inverse Fourier transform operation. Note 

in (25) that as a -*■ OO , the boundary condition becomes E x (t)=0.0, 

which is required for a perfect conductor. To determine Z ? '(t), 
the Laplace transform variable s=jw is used in (23) to obtain 
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Using the Laplace transform pair [11] 


1 



= sr 1 



(28) 


where the eg -1 denotes the inverse Laplace transform operation; and 
Z s ' (t) is then determined to be 


Z s (t) 






\J 7rat 

0, 


t >0 


t <0 


(29) 


This is the required time domain surface impedance function. 
Substituting (29) into (25) and discretizing space and time 
according to (15) gives 


EX n (k+l) 


\J 7ra(n£t) 


5(n6t) 


HY n (k+1/ 2 ) 


(30) . 


Substituting (30) into (16) yields 


-M 0 Sz 


d(nSt) 


HY n (k + 1/ 2 ) 


M 


\J 7ra(n<St) 


d (n£t) 


HY n (k+1/ 2 ) 


- EX n (k) 
(31) 


The convolution in (31) is expressed as a summation to obtain 


-H 0 6z 


d (n<$t) 




n-1 

HY n (k + 1/2 ) 

- 

M<St 


"N 

m=0 


d ( (n-m) St) 


(32) 


where Z 0 (m) is the discrete impulse response. The discrete impulse 
response is obtained by assuming the fields are piecewise constant 
in time as 


Z 0 (m) 


m+1/2 



(33) 
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If m=0 , the lower limit in (33) is 0. Approximating the time 
derivatives on the magnetic fields in (32) as finite differences 
results in 


[HY n+1/2 (k + l/2) - HY n1/2 (k + 1/ 2 ) ] = -Z 1 i:[HY n - ffl+1/2 (k + l/2) -HY n ' fl " V2 (k + 1/ 2 ) ] Z 0 (m) 

m=0 


+ -H- EX n (k) 




(34) 


where 


z i - 


m 0 5zn| 


ust 


7 TO 


(35) 


Equation (34) is suitable for computer implementation and includes 
the full convolution with all past field components. This full 
convolution would be impractical for large three dimensional 
problems; thus it is desirable to obtain a more efficient 
implementation. The development of a recursive implementation is 
the subject of the following section. 

V. Recursive Implementation 

Recently, Luebbers et. al. [3] extended the FDTD technique to 
dispersive dielectric materials using a time domain susceptibility 
function for polar dielectrics. In that paper, the time domain 
susceptibility function was a decaying exponential which permits 
the convolution summation to be recursively updated, thus avoiding 
the need for the complete time history of field components. Upon 
further examination of (34), it is clear that if Z 0 (m) can be 
approximated by a series of exponentials, then the SIBC can be 
efficiently evaluated using recursion. Figure 2 shows Z 0 (m) versus 
m, and it is clear that it can be approximated by a series of 
exponentials. Z 0 (m) is approximated as 

(36) 

1=1 

where N is the number of terms in the approximation. One of the 
most accurate methods for obtaining an exponential approximation to 
an exact function or to a data set is Prony's method [10]. Figure 
2 also shows the Prony approximation to Z 0 (m) with N=10 and it is 
clear that N=10 provides an adequate approximation. Thus, using 
(36) with N=10 in (34) gives 
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10 n-1 r 

£ 53 [HY n " m+1/2 (k+1/2) -HY nHn ' 1/2 (k + 1/2) 
EX n (k) 


HY n+1/2 (k + l/2) =HY n 1/2 (k + 1/2 ) - 


a m 

a ; e ' + 


1 + Zq ( 0 ) j = i m -i 

s t 


where 


/i 0 <Sz (l + Z^oCO) ) 


10 


z 0 (°) = E a < 


(37) 


(38) . 


i =1 


The convolution can now be recursively updated (see Appendix) to 
give 


10 


HY n+1 / 2 (k + 1/2) = HY n ~ 1/2 (k + 1/2 ) J3 (k + 1/2 ) 


<St 


H Q Sz (1 + Z^Zq(0) ) 


1 + Z^Zq( 0) 
EX n (k) 


(39) 


Note that only one past value of magnetic field is required to 
update the convolution summation. 

VI. One Dimensional Demonstration 

To demonstrate the constant and recursive FDTD SIBC, (20) and 
(39) were implemented in a one dimensional total field FDTD code 
for the geometry shown in Figure 1. The problem space size is 3 01 
cells, the impedance boundary is located at cell 3 00, and the 
electric field is sampled at cell 299. The maximum frequency of 
interest for each problem was 10 GHz. The incident electric field 
is a Gaussian pulse with maximum amplitude of 1000 V/m and has a 
total temporal width of 256 time steps. The frequency response of 
the incident pulse contains significant information to 12 GHz. Two 
computations were made with a— 2.0 S/ro and <7=20.0 S/m. The loss 
tangents at 10 GHz are 3.599, 35.99, respectively. The 

permittivity and permeability for the lossy dielectric were those 
of free space. The cell size and time step were 750 urn and 2.5 
psec, respectively. A tie point of 5.0 GHz was chosen for the FDTD 
constant SIBC. For each FDTD computation, a reflection coefficient 
versus frequency was obtained by first dividing the Fourier 
transform of the scattered field by the transform of the incident 
field at cell 299. The incident field was obtained by running the 
FDTD code with free space only and recording the electric field at 
cell 299. The scattered field is then obtained by subtracting the 
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time domain incident field from time domain total field. The 
results are compared with the analytic surface impedance reflection 
coefficient computed from 


R \ = ] Z S (fa)) - ^ol 

|Z s (a>) + n 0 


(40) 


where Z s (<i>) is given by (9) and T| 0 is the free space wave 
impedance. The phase of the FDTD reflection coefficient was 
corrected to account for the round trip phase shift of one cell 
since the FDTD reflection coefficient is computed from electric 
fields recorded one cell in front of the impedance boundary. 


The high conductivity surface impedance of (9) is an 
approximation to the general surface impedance for lossy 
dielectrics given by 


Z s (g>) 




'J (T + j(i)€ 


(41) . 


The advantage of using (9) over (41) for the FDTD SIBC 
implementation is that the resulting time domain impulse response 
is independent of the conductivity. The exponential approximation 
needs to be performed only once and not each time the conductivity 
is changed. 

Figures 3-4 show the FDTD constant and recursive SIBC 
reflection coefficient magnitude and phase results versus the 
analytic SIBC results for a=2.0 S/m. Notice the agreement between 
the curves is good, and the maximum error is about 0.02 at 10 GHz 
in Figure 3 . 

Figures 5-6 show the FDTD constant and recursive SIBC 
reflection coefficient magnitude and phase results versus the 
analytic SIBC results for a=20.0 S/m. Notice the agreement between 
the curves is excellent. 

Since the FDTD SIBC implementation is an approximation to an 
analytic SIBC, some amount of divergence between the SIBC curves 
and the analytic SIBC solution is to be expected with increasing 
frequency. As frequency increases, the effective number of cells 
per wavelength decreases and the FDTD SIBC becomes a rougher 
approximation to the analytic SIBC. To observe this error trend, 
the same one-dimensional test problems as above (using the 
dispersive SIBC only) were reevaluated with larger cell sizes equal 
to twice and four times the original cell size. This is equivalent 
to having 20 and 10 cells/A. 0 in the free space region, 
respectively. Figures 7 and 8 show the FDTD dispersive SIBC 


reflection coefficient magnitude and phase results versus the 
analytic SIBC for ct= 2.0 S/m using the original cell size and the 
larger cell sizes. Notice that for each increase in cell size, the 
agreement between the SIBC curve and the exact solution is reduced 
by a factor of two. This indicates that the error in the SIBC 
implementation is O(Sz) over the range of cell sizes examined here. 
The constant SIBC exhibited similar agreement reductions at the 5 
GHz tie point for larger cell sizes. 

VIII. Two Dimensional Demonstration 

As a practical application of the FDTD dispersive SIBC, 
frequency domain scattering width was computed from an infinite 
square cylinder for two scattering angles, 0=0.0 and 0=30.0 degrees 
using a full two dimensional TM scattered field code. The cylinder 
was 0.99 cm square and had parameters e=e 0 , M=M 0 f and cr=20.0 S/m. 
To illustrate the applicability of the SIBC, the cylinder was 
modeled in two ways. The first was a normal FDTD computation with 
a grid size of 10 cells/A (at 10 GHz) inside the conducting 
cylinder and the second was a SIBC computation with a grid size of 
10 cells/A. in free space (at 10 GHz) . Figure 7 shows the two 
dimensional field components and the cylinder dimensions (in cells) 
for the FDTD and SIBC computations. For the FDTD computation, the 
cylinder was modeled using 198 cells in the x and y directions, the 
cell size was 500 /urn, and the time step was 1.18 ps. For the SIBC 
computation, the cylinder was modeled using 32 cells in the x and 
y directions, the cell size was 0.003 m and the time step was 7.07 
ps. For both computations, a 100 cell border between the cylinder 
and the absorbing boundary was chosen, the total number of time 
steps was 1024, and an incident Gaussian pulse with total pulse 
width of 64 time steps was chosen. The near zone fields were 
transformed to far zone fields by a two-dimensional near zone to 
far zone transformation [12]. 

Figure 10 shows the scattering width magnitude versus 
frequency for a scattering angle of 0=0.0 degrees using the FDTD 
computation and the SIBC computation. Notice the good agreement 
over the entire frequency bandwidth for the dispersive SIBC. 

Figure 11 shows the scattering width magnitude versus 
frequency for a scattering angle of 0=30.0 degrees using the FDTD 
computation and the SIBC computation. Notice again the good 
agreement over the entire frequency bandwidth for the dispersive 
SIBC. 


IX. Summary 

One dimensional FDTD implementations of constant and 
dispersive surface impedance boundary conditions have been 
presented. The corresponding time domain impedance boundary 
conditions have been derived and their validity demonstrated by 
one-dimensional computation of the reflection coefficient at an 
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air-lossy dielectric interface at a single frequency and over a 
wide frequency bandwidth. The applicability of the SIBC to two- 
dimensional scattering problems was demonstrated by scattering 
width computation from an infinite square cylinder. For both the 
one and two dimensional cases, the dispersive FDTD results were 
shown to be in good agreement with exact results over the entire 
bandwidth. Considerable computational savings were illustrated and 
a recursive updating scheme was implemented which permits efficient 
application of a dispersive surface impedance boundary condition to 
practical scattering problems. 

Future extensions of this surface impedance concept currently 
under investigation are implementation in three dimensions, 
inclusion of surface curvature, dispersive dielectric and magnetic 
materials and thin material layers. 
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Appendix 

The purpose of this Appendix is to show how the discrete 
convolution of the SIBC in (37) can be done recursively. The 
convolution in (37) is 

10 n-1 r I 

££[ HY nHn+1/2 (k + l/2) - HY n " m 1/2 (k + 1/2 ) J a 1 e“ m (42) 

5=1 m=1 


Consider n=2, and (42) becomes 


^ ( H Y 3/2 (k + 1/2 ) - HY 1/2 (k + 1/2 ))a j e° 1 ’ 


i =1 


Now define 


(43) 


i}fj (k + 1/2) = (HY 3/2 (k + l/2) - HY 1/2 (k+1/2) ) a,.e“’ (44) 


Next for n=3, (42) becomes 

10 2 


YY (HY 7/2 ~ ra (k + 1/2) - HY 5/2 " m (k + l/2) )a,e 

l =1 m=1 


(45) 


Expanding (45) gives 


(hy 5/2 (k + 1/2) - HY 3/2 (k + 1/2 ))a,e*’ + 


i =1 

(. 


HY 3/2 (k + 1/2) - HY‘ ,c (k 


1 / 2 , 


:+l/2))a,.« 


2a 


(46) 


Substituting (44) into (46) we obtain 
10 


Y (hy 5/2 (k + 1/2) - HY 3/2 (k+l/2))a j e a ' + e“i|r- (k+1/2) 


i =1 


Equation (47) can be generalized for any time step n as 
1° n v 

Y \HY n ~ 1/2 (k + 1/2) - HY n 3/2 (k+1/2 ) ) a,e a ' + e 0 ’^ (k+1/2) 


(47) 


(48) 


with 
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i|r"(k + l/2) = (hY 0-172 (k + 1/2 ) -HY n ' 3/2 (k+l/2))a j e° t ' 

e a ’i|i n_1 (k + 1/2) 


(49) 


and 


i|r] (k + 1/2) = t?(k + l/2) = 0.0 


(50) 
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FIGURE TITLES 

1. Problem geometry showing one-dimensional FDTD grid and planar 
free space-conductor interface. 

2. FDTD dispersive SIBC discrete impulse response Z 0 (m) versus m 
and Prony approximation using 10 terms. 

3. Reflection coefficient magnitude versus frequency for normal 

incidence plane wave calculated for a=2 . 0 S/m using FDTD 

constant and dispersive SIBC and analytic solution. 

4. Reflection coefficient phase versus frequency for normal 

incidence plane wave calculated for a=2.0 S/m using FDTD 

constant and dispersive SIBC and analytic solution. 

5. Reflection coefficient magnitude versus frequency for normal 
incidence plane wave calculated for a=20.0 S/m using FDTD 
constant and dispersive SIBC and analytic solution. 

6. Reflection coefficient phase versus frequency for normal 

incidence plane wave calculated for a=20.0 S/m using FDTD 
constant and dispersive SIBC and analytic solution. 

7. Reflection coefficient magnitude versus frequency for normal 

incidence plane wave calculated for cr=2.0 S/m using FDTD 

dispersive SIBC with original and larger cell size and 

analytic solution. 

8. Reflection coefficient phase versus frequency for normal 

incidence plane wave calculated for a=2.0 S/m using FDTD 

dispersive SIBC with original and larger cell size and 

analytic solution. 

9. Two dimensional geometry for scattering width computations 
from an infinite square cylinder with a=20.0 S/m using normal 
FDTD and FDTD dispersive SIBC. 

10. Scattering width magnitude versus frequency at scattering 
angle 0=0.0 degrees from an infinite square cylinder with 
a=20.0 S/m using normal FDTD and FDTD dispersive SIBC. 

11. Scattering width magnitude versus frequency at scattering 
angle 0=30.0 degrees from an infinite square cylinder with 
a=20.0 S/m using normal FDTD and FDTD dispersive SIBC. 
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SUMMARY 

Surface impedance boundary conditions are employed 
to reduce the solution volume during the analysis of 
scattering from lossy dielectric objects. In a finite 
difference solution, they also can be utilized to avoid 
using small cells, made necessary by shorter wavelengths 
in conducting media, throughout the solution volume. 
This paper presents a one dimensional implementation for 
a surface impedance boundary condition for good 
conductors in the Finite Difference Time Domain (FDTD) 
technique. 

In order to illustrate the FDTD surface impedance 
boundary condition, we considered a planar air-lossy 
dielectric interface as shown in Figure 1. The incident 
field has polarization TE Z , and is propagating in the +z 
direction. The one-dimensional FDTD grid is also shown 
in Figure 1. To begin our implementation for a FDTD 
surface impedance boundary condition, we assume that the 
lossy dielectric has permittivity e , permeability /jl , and 
conductivity oi and that it is a good conductor. Thus, 
these constitutive parameters are real and satisfy the 
relation 


— >1 ( 1 ) 
toe 


where to is the radian frequency. We also assume that 
the radius of curvature is large compared to the maximum 
wavelength in the material and that the thickness of the 
material is large compared to the skin depth. 


The first order (or Leontovich) impedance boundary 
condition in the frequency domain is [1] 

E x {(x>) = T) g (u)Hy (<I>) (2). 

where is the surface impedance of the conductor. 

A superscript "t" is used in equation (2) to indicate the 
boundary condition is applied to the total field in free 
space. The frequency domain surface impedance for good 
conductors is 








(3) . 


Separating incident and scattered E x terms in (2) yields 
E*(a) =r\ g (a)Hy( (ji) -£*(«) (4). 

where the superscripts "s" and " i" are used to denote 
scattered and incident field components respectively. 
This equation is the required surface impedance boundary 
condition for scattered field components. The 

corresponding time domain expression involves a 
convolution integral and is given as 

E x {t) = ti -E x Ut) (5) 

where the ' * ' denotes convolution and the time domain 
surface impedance impulse response is given by 

r\' s (t) =F- 1 W a (u>)) = F-H ^4^ ) (6). 


We have approximated this time domain impulse response by 
a series of exponentials to obtain an efficient recursive 
updating scheme requiring only four running sum variables 
similar to [2]. Figures 2 and 3 show reflection 
coefficient comparison versus frequency for 
conductivities of 10.0 S/m and 50.0 S/m. The FDTD 
reflection coefficients are compared against the standard 
analytical solution. Note that the agreement is quite 
good for the entire frequency band. 

Overall, the surface impedance boundary condition 
implementation works well in eliminating the conductor 
volume from the solution space. This method has a 
distinct advantage over other possible implementations 
because the coefficients of the exponential approximation 
of the impulse response are independent of the 
conductivity of the scattering object and do not need to 
be reevaluated for different conductivities. Extensions 
of this surface impedance concept to two and three 
dimensions are currently under investigation. 
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Figure 1. Problem geometry showing incident plane wave, 
planar air-lossy dielectric interface and one 
dimensional FDTD grid. 



FDTD SURFACE IMPEDANCE COMPARISON 
10.0 S/m Conductivity 



Figure 2. Reflection coefficient magnitude versus 

frequency for normal incidence plane wave 
calculated for 10.0 S/m conductivity using 
FDTD surface impedance and analytical 
evaluation of high conductivity approximation. 
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Abstract 

Radar Cross Section (RCS) calculations for flat, perfectly conducting plates are 
readily available through the use of conventional frequency domain techniques such as 
the Method of Moments (MOM). However, if the plate is covered with a dielectric material 
that is relatively thick in comparison with the wavelength in the material, these frequency 
domain techniques become increasingly difficult to apply. In this paper, we present the 
application of the Finite Difference Time Domain (FDTD) technique to the problem of 
electromagnetic scattering and RCS calculations from a thin, perfectly conducting plate 
that is coated with a thick layer of lossless dielectric material. Both time domain and RCS 
calculations will be presented and discussed. 

I. Introduction 

The Finite Difference Time Domain (FDTD) technique has become increasingly 
popular in recent years for modeling electromagnetic scattering problems. It is based 
upon the time domain form of Maxwell’s equations, in which temporal and spatial 
derivatives are approximated by finite differences, and the electric and magnetic fields are 
interleaved spatially and temporally. Transient scattering behavior is easily examined and 
through the use of non-sinusoidal plane wave excitation, wideband frequency results can 
be obtained. The technique was first proposed by Yee [1] in 1966 and is inherently 
volumetric, which makes it ideal for modeling volumetric scatterers. Thin scatterers can 
easily be accommodated, and recently the technique has been expanded to include 
dispersive materials [2], plasmas [3], and chiral materials [4]. Through the use of a near 
to far zone transformation [5], far zone scattered fields (and thus RCS data) are readily 
available. This paper presents time domain scattering and RCS calculations over 0-3 GHz 
for several incidence angles from a thin, perfectly conducting (PEC) plate that is coated 
with a uniform lossless dielectric layer. 

II. Problem Description 

The scattering problem was a 3X by 6* (at 3 GHz) perfectly conducting plate that 
was coated with a 5 cm thick lossless dielectric layer with relative dielectric constant of 
€ r =4.0. Figure 1 shows the problem geometry with the dielectric layer on top of the plate. 

The wavelength (at 3 GHz) inside the dielectric layer is a. 0 /^7T =5 -° cm ’ where is the 

free space wavelength at 3 GHz. Thus, the dielectric coating is relatively thick at U. The 
spatial increment (cell size) was chosen to be 1 cm, which provides a spatial resolution 
of 5 cells/A. inside the dielectric coating and 10 cells/A. 0 in free space. 

The problem space size was chosen to be 61 by 121 by 49 cells in the x. y and 
z directions respectively. The plate was centered within the problem space in the x and 
y directions. The plate was positioned low in the problem space in the z direction 
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to allow any specular reflections multiple encounters with the outer radiation boundary 
condition (ORBC). 

The plate was constructed with 30 by 60 by 5 cells in the x, y and z directions for 
the dielectric coating, and with 30 by 60 by 0 cells for the PEC plate. The dielectric 
coating was constructed first, and the PEC plate was constructed on the bottom of the 
dielectric layer to avoid any air gaps within the scatterer. A 15 cell and 30 cell border on 
each side of the scatterer in the x and y directions provided adequate margin for the near 
to far zone transformation integration surface and for the ORBC. 

A 6 -polarized, Gaussian pulse incident plane wave with a maximum amplitude of 
1000 V/m and a total temporal width of 128 time steps was chosen. The time step was 

0. 0192 ns and the total number of time steps was 2048. 

III. Computations and Discussion 

Calculations were made at incidence angles 0 =0.0, 0 =60.0, 0 =85.0 and 0 =90.0 
degrees for both an uncoated and coated plate. The incidence angle was taken from the 
+z axis, the 0 incidence angle was <p =0.0 degrees for ail computations, and the far field 
computations were for backscatter only. The 0-polarized scattered and incident fields 
were then transformed to the frequency domain via an FFT and RCS was determined. 
Each computational problem required slightly more than one hour of CPU time on a Cray 
Y-MP supercomputer. 

Figure 2 shows the 0-polarized time domain far zone scattered electric field for 
incidence angle 0 =0.0 degrees. Note for the coated plate the early reflection from the 
top edge of the dielectric layer. Also note the time domain response for the coated plate 
is not markedly different from the uncoated plate except for some additional '‘ringing“ due 
to energy being confined within the dielectric coating. Figure 3 shows the RCS 
computations versus frequency again for both the uncoated and coated plate, and the 
RCS for both cases does not differ substantially. 

Figure 4 shows the 0-polarized time domain far zone scattered field for incidence 
angle 0 =60.0 degrees. Note the time responses differ more substantially for this case 
as more energy is being confined within surface wave modes of the dielectric layer. 
Figure 5 shows the corresponding RCS. Note the small peaks that have appeared in the 
RCS for the coated plate. We postulate these peaks correspond to surface wave modes 
that have been excited and radiate energy to the far field. As an approximation, we 
computed cutoff frequencies for waveguide modes for an infinite, dielectric covered 
ground plane according to Balanis [6]. These waveguide modes and corresponding 
cutoff frequencies are tabulated in Table 1. Examining the RCS of the coated plate, it is 
easily seen that a peak in the RCS is in close proximity to each cutoff frequency from 
Table 1. The peaks in the RCS are not located exactly at the cutoff frequencies of Table 

1, probably due to the finite size of the plate. 
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Mode 

Cutoff Frequency (GHz) 

TM 0 

0.00000 

TE, 

0.86589 

tm 2 

1.73718 

te 3 

2.59767 


Table 1 . Modes and cutoff frequencies for 
an infinite ground plane covered with a 
5 cm thick dielectric layer of e r =4.0. 

Figure 6 shows the 0 -polarized time domain far zone scattered field for incidence 
angle 0 =85.0 degrees. Also shown in Figure 6 is the far zone scattered field for a 5 cm 
thick 3A. by 6A dielectric layer. Since this incidence angle is near grazing, we expect to 
see little scattered field for the uncoated plate. Examining the time responses in Figure 
6, the uncoated plate response is indeed quite small in comparison to the dominant 
coated plate response. The dielectric layer response has the same general form as the 
coated plate response but is smaller in magnitude. Figure 7 shows the corresponding 
RCS. Note the large peaks and lobing structure for the coated plate and dielectric layer 
RCS. These can be attributed to radiation from surface wave modes of the dielectric 
cavity. 


Figure 8 shows the 0 -polarized time domain far zone scattered field for incidence 
angle© =90.0 degrees. The response for the zero-thickness uncoated plate is zero, while 
the coated plate time response does not differ substantially from that for incidence angle 
0 =85.0 degrees. Figure 9 shows the RCS for the coated plate only, and it is also similar 
to the coated plate RCS for incidence angle 0 =85.0 degrees. 

IV. Conclusions 

In this paper, the FDTD technique has been applied to model electromagnetic 
scattering from a perfectly conducting plate coated with a uniform, lossless dielectric 
layer. Time domain scattering results and frequency domain Radar Cross Section 
computations were presented and discussed. Large peaks in the RCS were founc for the 
coated plate at large incidence angles (near grazing) due to energy being radiated from 
surface wave modes of the dielectric layer. 
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The next step would be to provide a rigorous analytical treatment of the problem 
of a dielectric layer on a finite sized plate (ground plane) and to derive the surface wave 
structure of the layer and the far field scattering pattern. To the best of the authors’ 
knowledge, no such treatment has yet been presented. Results obtained from a rigorous 
theoretical treatment would then be used for comparison with the FDTD scattering 
computations and measured data. 
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Figure 1 . Problem geometry showing problem space size, plate size and plate position. 




THETA-POLARIZED BACKSCATTER 
3 BY 6 WAVELENGTH PLATE. THETA=0.0 



Figure 2. Monostatic, far zone, 0-polarized scattered field 
for uncoated and coated 3k by 6X plate at scattering angle 
6 =0.0 degrees. 


THETA- POLARIZED RADAR CROSS SECTION 
3 BY 6 WAVELENGTH PLATE, THETA =0.0 



Figure 3. Monstatic Radar Cross Section versus 
frequency for an uncoated and coated 3A. by 6A. plate at 
scattering angle 0 =0.0 degrees. 
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Figure 4. Monostatic, far zone, 0-polarized scattered field 
for uncoated and coated 3X by 6A. plate at scattering angle 
0 =60.0 degrees. 



Figure 5. Monostatic Radar Cross Section versus 
frequency for an uncoated and coated 3k by 6A. plate at 
scattering angle 0 =60.0 degrees. 
















Figure 6. Monostatic, far zone, 0 -polarized scattered field 
for uncoated and coated 3k by 6A. plate and for 5 cm thick 
3A. by 6A. dielectric layer at scattering angle 0=85.0 
degrees. 



Figure 7. Monostatic Radar Cross Section versus 
frequency for an uncoated and coated 3k by 6* plate and 
for a 5 cm thick 3k by 6k dielectric layer at scattering 
angle 0 =85.0 degrees. 
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THETA-POLARIZED BACKSCATTER 
3 BY 6 WAVELENGTH PLATE, THETA=90.0 
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Figure 8. Monostatic, far zone, 0 -polarized scattered field 
for uncoated and coated 3A by 6 k plate at scattering angle 
0=90.0 degrees. 
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Fgure 9, Monostatic Radar Cross Section versus 
frequency for a coated 3A. by 6X plate at scattering angle 
0 =90.0 degrees. 
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Abstract 

Radar Cross Section (RCS) calculations for flat, perfectly conducting plates are 
readily available through the use of conventional frequency domain techniques such as 
the Method of Moments. However, if time domain scattering or wideband frequency 
domain results are desired, then the Finite Difference Time Domain (FDTD) technique is 
a suitable choice. In this paper, we present the application of the Finite Difference Time 
Domain (FDTD) technique to the problem of electromagnetic scattering and RCS 
calculations from a thin, perfectly conducting plate for a conical cut in the scattering angle 

0 . RCS calculations versus angle 0 will be presented and discussed. 

1. Introduction 

The Finite Difference Time Domain (FDTD) technique has become increasingly 
popular in recent years for modeling electromagnetic scattering problems. It is based 
upon the time domain form of Maxwell’s equations, in which temporal and spatial 
derivatives are approximated by finite differences, and the electric and magnetic fields are 
interleaved spatially and temporally. Transient scattering behavior is easily examined and 
through the use of non-sinusoidal plane wave excitation, wideband frequency results can 
be obtained. The technique was first proposed by Yee [ 1 ] in 1966 and is inherently 
volumetric, which makes it ideal for modeling volumetric scatterers. Thin scatterers can 
easily be accommodated, and recently the technique has been expanded to include 
dispersive materials [2], plasmas [3], and chiral materials [4]. Through the use of a near 
to far zone transformation [5], far zone scattered fields (and thus RCS data) are readily 
available. This paper presents RCS calculations at several frequencies of a conical cut 
in scattering angle 0 for a 3.5A. by 2X perfectly conducting plate. Results obtained with 
FDTD computations are compared with results obtained with the ESP4 electromagnetic 
code. 

II. Problem Description 

This particular scattering problem first came to the authors’ attention when Dr. Woo 
indicated he was obtaining significant discrepancies between measurements and 
numerical analyses using several Method of Moments electromagnetic codes. He then 
suggested that a scattering analysis by the Penn State FDTD code may provide some 
insight as to where problems may exist. 

The scattering problem was a 3.5k by 2k (at 3 GHz) perfectly conducting plate that 
was oriented within the FDTD solution space as shown in Figure 1. The free space 
wavelength (at 3 GHz) is k 0 = 10.0 cm. The spatial increment (cell size) was chosen to 
be 1 cm, which provides a spatial resolution of 10 cells/X 0 in free space. 
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The problem space size was chosen to be 66 by 41 by 49 cells in the x, y and z 
directions respectively. The plate was centered within the problem space in the x and y 
directions. The plate was positioned low in the problem space in the z direction to allow 
any specular reflections multiple encounters with the outer radiation boundary condition 
(ORBC). 

The plate was constructed with 35 by 20 by 0 cells in the x, y and z directions for 
the PEC plate. Thus the physical plate size was 35 cm by 20 cm. A 15 cell and 10 cell 
border on each side of the scatterer in the x and y directions provided adequate margin 
for the near to far zone transformation integration surface and for the ORBC. 

A 0 -polarized, Gaussian pulse incident plane wave with a maximum amplitude of 
1000 V/m and a total temporal width of 128 time steps was chosen. The time step was 
0.0192 ns and the total number of time steps was 1024. 

III. Computations and Discussion 

At the recommendation of Dr. Woo, calculations were made for 0 =80.0 degrees; 
and the angle <p was varied from 0.0 to 10.0 degrees in steps of 0.5 degrees. The 
incidence angles 0 and <p were taken from the +z and +x axes respectively and the 
incident field was <p -polarized for all computations. The far field computations were for 
backscatter only. For each incidence angle, the 0 -polarized (co-pol) and 0 -polarized 
(cross-pol) scattered and incident fields were then transformed to the frequency domain 
via an FFT and RCS was determined. From each RCS data file, the 0 -polarized and 0 - 
polarized RCS at 3 GHz for each angle 0 were chosen and then entered into separate 
data files of RCS versus angle 0 . Each incidence angle computation required 1.5 hours 
on a 486/25 personal computer. 

Figure 2 shows a sample time domain co-pol backscatter for 0 =80.0 degrees and 
0 =5.0 degrees. Figure 3 shows a sample time domain cross-pol backscatter for 0 =80.0 
degrees and <p =5.0 degrees. 

Figure 4 shows the co-pol radar cross section versus scattering angle 0 . Note the 
extreme disagreement between the FDTD result and the ESP4 result. 

Figure 5 shows the cross-pol radar cross section versus scattering angle 0 . Note 
the change in scale of the cross section and the relatively good agreement between the 
FDTD result and the ESP4 result. Since the co-pol RCS results for FDTD and ESP4 were 
very much different and the cross-pol RCS results were similar, this indicates a possible 
problem with either the FDTD code or the ESP4 code. 
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IV. Investigation 

The original calculations for ESP4 were made at 3 GHz with a plate size of 35 cm 
by 20 cm. These results were compared with FDTD RCS frequency domain data at 3 
GHz. In order to gain some insight as to possible problems with the FDTD or ESP4 
electromagnetic codes, three additional RCS computations were made with ESP4 at 
frequencies of 2.25 GHz, 2.47 GHz and 2.70 GHz with 0 =80.0 degrees and a#-polarized 
incident field. These frequencies were chosen as approximate frequencies for two peaks 
(2.25 and 2.70 GHz) and one null (2.47 GHz) in the co-pol frequency domain RCS from 
Figure 6. The data points from the resulting FDTD RCS files were chosen for each 
frequency and entered into a separate data file for comparison with ESP4 results. 

Figure 6 shows the co-pol radar cross section versus frequency and incidence 
angle #. Figure 7 shows the cross-pol radar cross section versus frequency and 
incidence angle#. 

Figures 8 and 9 show the co-pol and cross-pol radar cross section versus 
scattering angle # for both FDTD and ESP4 at 2.25 GHz. Note also the extreme 
disagreement between the co-polarized RCS results and the good agreement between 
the cross-polarized RCS results. 

Figures 10 and 11 show the co-pol and cross-pol radar cross section versus 
scattering angle # for both FDTD and ESP4 at 2.47 GHz. Again, note the disagreement 
between the FDTD and ESP4 co-polarized RCS and the good agreement between the 
cross-polarized RCS results. 

Figures 12 and 13 show the co-pol and cross-pol radar cross section versus 
scattering angle # for both FDTD and ESP4 at 2.7 GHz. Note the co-polarized RCS 
results are in slightly better agreement and the cross-polarized results again are in good 
agreement. 

V. Conclusions 

In this paper, the FDTD technique has been applied to model electromagnetic 
scattering in a conical cut from a perfectly conducting plate. Radar Cross Section versus 
scattering angle# was presented and discussed. The# -polarized FDTD and ESP4 RCS 
results were very dissimilar, while the 0 -polarized RCS results were in good agreement. 
Thus a problem may exist with one or both of the FDTD or ESP4 electromagnetic codes. 
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Figure 1. Problem geometry showing FDTD solution space size, and the size 
and spatial placement of the perfectly conducting plate. Dimensions are in cells, 
with each cell being a 1 cm cube. 



Figure 2. Sample time domain far zone, backscattered, 
co-polarized electric field from a 35 cm by 20 cm 
perfectly conducting plate. The scattering angles are 
9 =80.0 , 0 =5.0 degrees. 



Figure 3. Sample time domain far zone, backscattered, 
cross-polarized electric field for a 35 cm by 20 cm 
perfectly conducting plate. The scattering angle is 
0 =80.0 , 0 =5.0 degrees. 
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Figure 4. Co-polarized radar cross section versus angle 0 for 
a 35 cm by 20 cm perfectly conducting plate at 3.00 GHz 
using FDTD and ESP4. 
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Figure 5. Cross-polarized radar cross section versus angle <p 
for a 35 cm by 20 cm perfectly conducting plate at 3.00 GHz 
using FDTD and ESP4. 
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Figure 7. Cross-polarized radar cross section versus frequency and incidence angle <t> from a 35 cm by 20 cm perfectly 
conducting plate using FDTD. 






Figure 8. Co-polarized radar cross section versus angle 0 
from a 35 cm by 20 cm perfectly conducting plate at 2.25 GHz 
using FDTD and ESP4. 
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Figure 9. Cross-polarized radar cross section versus angle 0 
from a 35 cm by 20 cm perfectly conducting plate at 2.25 GHz 
using FDTD and ESP4. 





Figure 10. Co-polarized radar cross section versus angle 0 
from a 35 cm by 20 cm perfectly conducting plate at 2.47 GHz 
using FDTD and ESP4. 
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Figure 1 1 . Cross-polarized radar cross section versus angle 
<p from a 35 cm by 20 cm perfectly conducting plate at 2.47 
GHz using FDTD and ESP4. 
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Figure 12 . Co-polarized radar cross section versus angle 0 
from a 35 cm by 20 cm at 2.70 GHz using FDTD and ESP4. 



Figure 13. Cross-polarized radar cross section versus angle 
0 from a 35 cm by 20 cm perfectly conducting plate at 2.70 
GHz using FDTD and ESP4. 
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Abstract 

Thin sheets of resistive or dielectric material are commonly 
encountered in radar cross section calculations. Analysis of 
such sheets is simplified by using sheet impedances. In this 
paper it is shown that sheet impedances can be modeled easily and 
accurately using Finite Difference Time Domain (FDTD) methods. 
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Introduction 


In [1] a review of various approximate boundary conditions is 
given, including several for thin sheets and layers. These are 
applicable to sheets which are thin relative to the free space 
wavelength, so that they can be approximated by an electric 
current sheet. If the thin sheet is primarily conductive the 
sheet impedance will be resistive, as is the case for resistance 
cards. A thin lossless dielectric sheet will have a purely 
reactive sheet impedance, while in general the sheet impedance 
will be complex. These sheets are characterized by a 
discontinuity in the tangential magnetic field on either side of 
the sheet but no discontinuity in tangential electric field. 

This continuity, or single valued behavior of the electric field, 
allows the sheet current to be expressed in terms of an impedance 
multiplying this electric field. 

Approach 

The sheet impedance can be defined in several ways. A 
convenient definition can be obtained by combining eqs. (3.3) and 
(3.5) of [2] 


= oT + joe 0 (e r -l) T (i) 

with 

= l/r 3 (2) 

where Y s is the sheet admittance, Z s the sheet impedance, a and e r 
the conductivity and relative permittivity of the sheet material, 
T the sheet thickness, and e Q the free space permittivity. 

Let us now consider how to incorporate this approximation 
into the FDTD method. The surface impedance approximation 
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requires the impedance sheet to be small compared with the free 
space wavelength. In most FDTD calculations the FDTD cell size 
(Yee [3] cells are used here) must be on the order of l/io 
wavelength or less for reasonably accurate results. Scattering 
from an infinitesimally thin perfectly conducting plate was has 
been calculated by approximating the plate as being one FDTD cell 
thick with good results [4]. If it is assumed that the same 
approach can be applied to infinitesimally thin impedance sheets, 
then the plate thickness T in (l) merely becomes the thickness of 
the FDTD cell, and the conductivity and/or relative permittivity 
to be used in the FDTD calculations are merely adjusted in 
accordance with (1) to give the desired sheet impedance. Note 
that the FDTD cell dimension need not correspond to the thickness 
of the actual physical plate. The FDTD cell thickness is used 
only to determine the conductivity and relative permittivity of 
the FDTD electric field location so that the desired sheet 
impedance is approximated. Note also that, even if the 
wavelength in the material forming the impedance sheet is much 
smaller than a free space wavelength, the FDTD cell size need not 
be correspondingly reduced. 

Demonstration 


The first demonstration will consist of calculating the far 
zone backscatter from a 29 x 29 cm flat plate of sheet impedance 
Z s = 500 n. The FDTD calculations will use cubical Yee cells 
with 1 cm edges. Using T = l cm, the corresponding FDTD 

conductivity is a = 0.2 S/m. The FDTD calculations shown in 

Figures 1-8 are all made with the plate modeled by setting the 
conductivity to 0.2 S/m for x and y polarized electric field 

locations corresponding to single z dimension index over a range 

of x and y dimension indices to model the plate. The FDTD 
approach used and the transformation to the far zone is described 
in [4]. The problem space size, orientation and position of the 
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plate, incident Gaussian pulse plane wave, and time step size are 
also consistent with those in [4]. 

Figure 1 shows the far zone backscattered electric field for 
a Gaussian pulsed plane wave normally incident on the plate. In 
Figure 2 this result is Fourier transformed, converted to cross 
section, and compared with results using the Method of Moments 
[2]. The agreement is guite good, with the approximately 20 dB 
reduction in radar cross section relative to a perfectly 
conducting plate of the same size [4] consistently predicted by 
both methods. 

In Figures 3-8 the same plate geometry and composition is 
considered but for non-normal incidence. The plate is 
perpendicular to the z axis, with edges parallel to the x and y 
axes, and the plane wave is incident from 0=45, 0=30 degrees. 
Figures 3-5 show the co-polarized backscatter far zone electric 
field for <p and 0 polarizations and the cross-polarized 
backscatter as well. in Figures 6-8 these time domain results 
are Fourier transformed and converted to radar cross section for 
comparison with Moment Method [2] results. Again the agreement 
is quite good, except at the highest frequencies considered. 

These results indicate that perhaps 12 cells/ wavelength are 
required for good accuracy for off-normal incidence. Comparing 
the results in Figure 6 with those in Figure 5 of [4], it is 
clear that changing from a perfectly to a finitely conducting 
piste changes the scattering level and frequency behavior, and 
that the FDTD and Moment Method results agree quite well on these 
effects. 

In Figure 9 both FDTD and Moment Method [4] results for 
scattering by a plate with a complex sheet impedance are shown. 
The sheet impedance is determined by applying eqs. (1,2) with 
conductivity 0.25 S/m, relative permittivity 3.0, and thickness l 
cm. , corresponding to the FDTD parameters used. Again the plane 
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wave is a Gaussian pulse incident from 0=45, 0=30 degrees. The 
FDTD results agree with the Moment Results for frequencies up to 
about 12 cells/wavelength. 

The final result is for a plate with edge treatment. For 
this demonstration a 21 x 21 cm thin perfectly conducting plate 
is given a 4 cm border of sheet impedance Z s = 500 n, resulting 
in a square plate 29 x 29 cm. This edged plate is modeled in 
FDTD by setting x and y polarized electric field locations for a 
single z dimension index as being either perfect conductor for 
the central portion of the plate or with a conductivity of 0.2 
S/m for the edges. The ESP4 calculations were made with a 
central perfectly conducting plate surrounded by 4 plates of 
sheet impedance Z s = 500 fl attached to the central plate using 
overlap modes. The results are compared in Figure 10 with 
excellent agreement between the two methods, both showing a 
significant difference due to the edge treatment when compared 
with the results of Figure 6 . 


Conclusions 


The ability of the FDTD method to easily and accurately 
model scattering by sheet impedances was demonstrated by 
comparing FDTD results for scattering from flat plates modeled 
using sheet impedances with Method of Moment results. The 
approach described here is directly applicable to the Yee cell, 
and demonstrated good accuracy for frequencies up to 
approximately 12 cells per wavelength. 
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Figure Titles 


1. Co-Polarized far zone electric field vs time scattered by a 
29 x 29 cm flat plate of sheet impedance 500 ohms for a 0- 
polarized normally incident Gaussian pulse plane wave 
computed using FDTD. 


2 . 


3. 


4. 


Radar cross section for a 29 x 29 cm flat plate of sheet 
impedance 500 ohms, normal incidence, obtained from FDTD 
results of Figure 1 compared with Moment Method [2] results. 


Co-Polarized far zone electric field vs time scattered by a 
29 x 29 cm flat plate of sheet impedance 500 ohms for a 0- 
polarized incident Gaussian pulse plane wave from 0=45, 0=30 
degrees computed using FDTD. 


Co-Polarized far zone electric field vs time 
29 x 29 cm flat plate of sheet impedance 500 
polarized incident Gaussian pulse plane wave 
degrees computed using FDTD. 


scattered by a 
ohms for a 0- 
from 0=45, 0=30 


5. Cross-Polarized far zone electric field vs time scattered by 
a 29 x 29 cm flat plate of sheet impedance 500 ohms for a <p- 
polarized incident Gaussian pulse plane wave from 0=45, 0=30 
degrees computed using FDTD. 


6. Co-Polarized radar cross section for a 29 x 29 cm flat plate 
of sheet impedance 500 ohms, 0=45, 0=30 degree incidence, 0- 
polarized, obtained from FDTD results of Figure 3 compared 
with Moment Method [2] results. 

7. Co-Polarized radar cross section for a 29 x 29 cm flat plate 
of sheet impedance 500 ohms, 0=45, 0=30 degree incidence, 0- 
polarized, obtained from FDTD results of Figure 4 compared 
with Moment Method [2] results. 

8. Cross-Polarized radar cross section for a 29 x 29 cm flat 
plate of sheet impedance 500 ohms, 0=45, 0=30 degree 
incidence, obtained from FDTD results of Figure 5 compared 
with Moment Method [2] results. 


9. Co-Polarized radar cross section for a 29 x 29 cm flat plate 
of sheet impedance corresponding to conductivity of 0.25, 
relative permittivity of 3.0, and thickness 1 cm, for 0=45, 
0=30 degree 0-polarized incident plane wave calculated using 
FDTD and compared with Method of Moments [ 2 ] . 


10. Co-Polarized radar cross section for a 21 x 21 cm perfectly 
conducting flat plate with a 4 cm 500 ohm edge treatment on 
all sides (total plate size 29 x 29 cm) for 0=45, 0=30 
degree 0-polarized incident plane wave calculated using FDTD 
and compared with Method of Moments [ 2 ] . 
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Abstract 

Determining transient electromagnetic fields in antennas with 
nonlinear loads is a challenging problem. Typical methods used 
involve calculating frequency domain parameters at a large number 
of different frequencies, then applying Fourier transform methods 
plus nonlinear equation solution techniques. If the antenna is 
simple enough so that the open circuit time domain voltage can be 
determined independently of the effects of the nonlinear load on 
the antenna current (an infinitesimal dipole, for example) , time 
stepping methods can be applied in a straightforward way. In 
this paper transient fields for antennas with more general 
geometries are calculated directly using Finite Difference Time 
omam methods. In each FDTD cell which contains a nonlinear 
load, a nonlinear equation is solved at each time step. As a 
test case the transient current in a long dipole antenna with a 
nonlinear load excited by a pulsed plane wave is computed using 
this approach. The results agree well with both calculated and 
measured results previously published. The approach given here 
extends the applicability of the FDTD method to problems 
involving scattering from targets including nonlinear loads and 
materials, and to coupling between antennas containing nonlinear 

loads. It may also be extended to propagation through nonlinear 
materials. 
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INTRODUCTION 


Calculating the transient electromagnetic fields in antennas 
and scatterers containing nonlinear loads or material is a 
difficult problem. The traditional approach has been to apply 
frequency domain methods at a large number of harmonic 
frequencies. Sarkar and Weiner [1] used this approach in 
combination with a Volterra series analysis to determine 
scattering from antennas with nonlinear loads. Liu and Tesche 
[2,3] separated the problem into linear and nonlinear portions, 
calculated wideband frequency domain characteristics for the 
linear portion of the problem, transformed this to the time 
domain, and then solved the nonlinear portion by time marching. 
They achieved good agreement with measurements using this 
approach, and their results presented in [3] will be used to 
validate the results obtained in this paper. However, obtaining 
the frequency domain results for a complicated antenna at the 
large number of frequencies involved may be quite tedious and 
involve significant computer resources. And this approach cannot 
be extended to situations involving scattering from bulk regions 
of nonlinear materials, or to propagation through nonlinear 
materials as can the FDTD approach presented here. 

For the simpler situation where the time domain open circuit 
voltage on the antenna terminals can be determined independently 
from the effects of the nonlinear load the frequency domain 
portion of the approach of Liu and Tesche can be dispensed with, 
and the open circuit voltage can be used as the input to a 
nonlinear circuit model of the load, with the resulting problem 
solved directly using time marching or other methods applied to 
nonlinear circuits. This approach was used by Kanda [4], who 
also gives an excellent review of previous work in this area. 

Finally, Schuman [5] applied a time domain method of moments 
approach to a thin straight wire. However, his method appears 
difficult to apply to more general geometries. 

In this paper the Finite Difference Time Domain (FDTD) 
approach will be extended to include nonlinear lumped loads. The 
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reader is assumed to have some familiarity with this method. The 
literature is extensive, with some representative papers included 
in the following references [6-8], The authors are using the 
scattered field formulation of [8], but with linear time 
differencing. 

APPROACH 

To illustrate the method, let us consider the specific 
example of interest, taken from [3]. A wire dipole with half 
length 0.6 meters and a diameter of 0.81 mm is loaded at its 
midpoint, as shown in Figure 1. The dipole is located parallel 
to the z axis. The FDTD cell at the center of the wire is used to 
model the lumped load. As described in [3], this load is two 
diodes in series with a 100 ohm resistor (the actual measurements 
were made using a single diode at the base of a monopole) . The 
total diode junction capacitance of 0.5 pF must also be included 
in the model for accurate results at the frequencies contained in 
the pulse. 

In order to describe the approach used to model the diode 
circuit in FDTD, let us first consider an approach to 
approximating a linear lumped load consisting of a capacitor in 
parallel with a resistor (conductance) in an FDTD cell. Starting 
with 


V X H = e 


dE 

at 


+ oE 


(1) 


where H, E, e and a are the magnetic and electric fields, the 
permittivity, and the conductivity, and following the Yee [6] 
approach for discretizing space and time, with t=nAt, x=Iax, 
y=jAy, z=Kaz, eq. (1) becomes, for the z component of electric 
field in a particular FDTD cell. 
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(2) 


(V x H n+1/2 ) z 



a 



where e and a pertain to the particular cell location of E z and 
the superscripts denote the time reference of the particular 
component. Usually in applying FDTD this equation is solved for 
e"* 1 in terms of the previous time values of E n and H 1 **, with the 
curl H term computed using spatial finite differences. Instead, 
consider the physical meanings of the terms. The curl H term 
gives the total current density flowing in the cell surrounding 
the electric field component. The next term involving e and the 
time derivative of E is the displacement current density flowing 
through the cell in the z direction. The crE"* 1 term is the 
conduction current density flowing through the cell in the z 
direction. Using the cell dimensions ax, Ay and az, and assuming 
fields are constant across the cell, we can rewrite the above 
equation in terms of lumped elements, voltages, and currents as 

n+1 n 

AxAy(V x H n+1/2 ) = C Az — - + G Az e" +1 ( 3 ) 

z <St 


where now the first term is the total current flowing through the 
cell, C is the lumped "parallel plate" capacitance of the cell, 
and G is the lumped conductance in parallel with the capacitance. 
Note that one can identify az e"* 1 as the voltage across the cell. 
Clearly a lumped capacitance C can be equivalent to setting an 
appropriate value of the e of the cell based on the cell 
dimensions, and similarly for a lumped resistance and the a of 
the cell. Thus a lumped load that is a parallel combination of a 
capacitor and a conductance (resistance) can be modeled simply by 
sett ing the cell values of e and a appropriately. Equations (2) 
and (3) are interchangeable in terms of solving for e"* 1 . 
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However, one warning is that the cell permittivity cannot be 
set too low. FDTD in the form presented in this paper cannot in 
general model materials with an epsilon that is too small (much 
less than free space) without becoming unstable. Such materials 
can be modeled using FDTD modified for frequency dependent 
ma terials t 9 ], but with additional computational effort. If the 
conductance G (conductivity a) is great enough so that the 
displacement current term in eq. (2) (or (3)) can be neglected 
then this term may be dropped. Indeed, if (2) or (3) is solved 
for E and G (or equivalently a) is allowed to go to infinity, 
the correct result of E =0 is obtained. However, making G (or 
cr) small enough so that the conduction current is smaller than 
the displacement current through the cell filled with free space, 
but nevertheless neglecting this displacement current term, will 
result in instabilities. A physical argument for this is that 
the capacitance of an FDTD cell cannot be made lower than the 
capacitance of the cell filled with free space by adding lumped 
elements in parallel with the cell capacitance. 

Now let us proceed to extend this approach to the circuit of 
interest, shown in Figure 2. In this figure the resistance R 
models the input resistance of the oscilloscope used to measure 
the current [3]. The capacitance C represents the capacitance of 
the free space FDTD cell plus the junction capacitance of the 
diode. The diode junction capacitance of the diode is not 
actually in parallel with the resistor, but this approximation 
simplifies the following derivation and does not appreciably 
affect the results for the circuit element sizes under 
consideration. Letting i d represent the current through the 
diode, which is also the total conduction current through the 
cell, we can solve (3) for i d obtaining 

i d = AxAy(VxH n+1/2 ) z - — (eJ +1 -EJ (4) 

At 
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Once i d is determined from (4), the diode voltage v d can be 
obtained from the equation 

2 v d + i d R = v c = Az e" +1 (5) 

Since the diodes are in the circuit, i d must also satisfy the 
nonlinear equations 


i d = l.OxlO -8 v d , v d <0 (6a) 

i d = 2 . 9x10 7 [exp(15v d ) -1], v d >0 (6b) 

where (6) are the nonlinear diode equations given in [3]. Since 
i d as given in both eqs. (4) and (6a, b) must be equal, a Newton- 
Raphson iteration method can be applied to solve for the e"* 1 
value which produces required equality. This was the approach 
taken in this paper. The convergence was very fast since an 
initial guess for E 1 ^ 1 of the previous value of E, E n , provided 
the Newton-Raphson iteration with a good starting value. 

DEMONSTRATION 

The dipole considered in [3] is excited by a pulsed plane 
wave. As shown in [3], this plane wave has a peak electric field 
strength of approximately 390 volts. For simplicity, this pulse 
has been approximated for our calculations by a Gaussian pulse. 
The FDTD excitation pulse used for the calculations in this paper 
is shown in Figure 3. The actual pulse shown in [3] rings at a 
low amplitude out to about 3 ns, but this was neglected in the 
FDTD calculations. 

In order to provide the necessary temporal resolution FDTD 
cells were chosen as 0.006 m cubes, providing a time step of 
11.55 ps. For the FDTD Gaussian pulse used in this demonstration 
this allowed 64 time steps between the 1% of peak amplitude 
values. The transient currents given in [3] have a duration of 
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approximately 14 ns. The corresponding FDTD calculation was 1300 
time steps (allowing for the incident pulse to reach the dipole) , 
reguiring 200 minutes on a 25 MHz 486 PC clone running Lahey 
fortran. 

To approximate the wire dipole 200 FDTD cells were used. 
Since the wire diameter is smaller than the FDTD cell width, 
subcell modeling was used to adjust for this [10]. The FDTD 
problem space was 39 x 39 x 240 cells, and was terminated in 
second order Mur [11] absorbing boundaries. 

Three calculations of the current through the dipole load 
were made and compared with results calculated by Liu and Tesche 
[3]. For all FDTD results the total current flowing through the 
FDTD cell, determined by evaluating the curl of H around the cell 
containing the lumped load, is plotted. In the first calculation 
only the 100 ohm resistor (in parallel with the free space 
capacitance of the FDTD cell) was included, with the results 
shown in Fig. 4. The agreement with the results of Liu is quite 
reasonable. 

Next the diodes were added in series with the 100 ohm 
resistor, and the FDTD cell capacitance C, shown in Fig 2, was 
set to 0.5 pF to model the combined diode junction capacitance. 
FDTD results for the two cases considered in [3] were then 
calculated. These are shown in Figs. 5 and 6, and differ only in 
that the pulse initially forward biases the diode for the results 
in Fig. 5, but reverse biases the diode in Fig. 6. Some points 
taken from calculated results in [3] are included in Figs. 5 and 
6 for comparison, but the interested reader on referring to [3] 
directly will see that the agreement between the FDTD results and 
both the calculated and measured results in [3] is excellent 
considering the different assumptions and approximations made in 
the analysis. 
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CONCLUSIONS 


In this paper an approach to model nonlinear lumped elements 
in the context of FDTD was presented. It was used to compute the 
transient current in a diode loaded long dipole antenna, and 
excellent agreement with previously published results was 
obtained. This capability extends the applicability of the FDTD 
method to a wide range of problems, including scattering from 
nonlinear loaded antennas and harmonic product generation by 
nonlinear elements. The method can be further extended to 
considering regions of nonlinear material, since the (nonlinear) 
material content of each FDTD cell can be specified 
independently. Thus extensions to scattering from targets 
containing nonlinear material, or to propagation through 
nonlinear media should be straightforward. In combination with 
FDTD methods for frequency dependent materials, it may also be 
possible to extend the method to model electromagnetic 
propagation through dispersive nonlinear materials, including 
soliton propagation through such media. 
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FIGURE TITLES 

Fig l. Pulsed plane wave incident on vertical wire dipole with 
nonlinear load. 

Fig 2. Approximate eguivalent circuit of lumped nonlinear 
load. Capacitance C includes both the FDTD cell 
capacitance and the diode junction capacitances. 
Resistance R simulates the input resistance of the 
oscilloscope used to measure the current [3]. 

Fig. 3 Gaussian pulse used in the FDTD calculations to 
simulate the pulse used by Liu et al in [3], 

Fig. 4 Transient current flowing through 100 ohm load and 

parallel FDTD cell capacitance at the terminals of the 
wire dipole calculated using FDTD and compared with 
calculated results of Liu et al [3]. 

Fig. 5 Total transient current flowing through equivalent 

circuit of Fig. 2 located at the terminals of the wire 
dipole calculated using FDTD and compared with 
calculated results of Liu et al [3]. Diode is 
initially forward conducting. 

Fig. 6 Total transient current flowing through equivalent 

circuit of Fig. 2 located at the terminals of the wire 
dipole calculated using FDTD and compared with 
calculated results of Liu et al [3], Diode is 
initially reverse conducting. 
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Abstract 


The Finite Difference Time Domain (FDTD) technique has been 
applied to a wide variety of electromagnetic analysis problems, 
including shielding and scattering. However, the method has not 
been extensively applied to antennas. In this short paper 
calculations of self and mutual admittances between wire antennas 
are made using FDTD and compared with results obtained using the 
Method of Moments. The agreement is quite good, indicating the 
possibilities for FDTD application to antenna impedance and 
coupling. 


I Introduction 

The Finite Difference Time Domain (FDTD) technique has had only 
limited application to antennas. This is somewhat surprising, 
since the geometrical and material generality of the method 
suggests that it might have significant application to antenna 
analysis, especially in situations where other structures, 
especially electromagnetically penetrable ones, are nearby. This 
is due to the relative ease with which the FDTD method 
accommodates modeling of volumetric electromagnetic interactions 
with materials as compared to the Method of Moments. 

Earlier work [1] has shown that the FDTD method could compute the 
self impedance of a wire antenna in three dimensions, however, 
the approach used in [1], plane wave incidence, did not lend 
itself to mutual coupling calculations. Accurate self -admittance 
FDTD results for two-dimensional antenna geometries were 
presented in [2]. In this paper we demonstrate both self and 
mutual admittance FDTD calculations for three dimensional wire 
antennas . 
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II Approach 


The test problem geometry is shown in Figure 1. Two wire dipoles 
of length 57 and 43 cm are parallel and separated by 10.5 cm. 

Both are center fed, and are symmetrically positioned. The goal 
is to determine the self admittance of the driven dipole and the 
mutual admittance between the two dipoles. 

The FDTD computations were made using a three dimensional 
computer code based on the Yee [3] cell, with second order Mur 
[4] absorbing boundaries. The problem space was chosen as 
61x51x80 cells, with the cell dimensions Ax=Ay=0.5 cm, Az=1.0 cm. 
Making the two transverse dimensions smaller results in a greater 
length to diameter ratio, so that a thin wire Moment Method code 
may be used to provide comparison results over a wider band of 
frequencies. Thinner wires may be modeled in FDTD using sub-cell 
methods [5,6] . 

For the FDTD calculations the longer dipole is fed at the center 
with a Gaussian pulse of 100 volts maximum amplitude that reached 
its 1/e amplitude in 16 time steps. The time steps were 11.11 
picoseconds, the Courant stability limit for the cell size 
chosen. 

During the progress of the FDTD calculations the currents at the 
center of each dipole were saved for each time step. They were 
computed by evaluating the line integral of the magnetic field 
around the dipoles at the center. Along with the applied 
Gaussian voltage pulse the currents were Fourier transformed to 
the frequency domain. Then, based on the admittance parameter 
equations 


h = v i 
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and taking V 1 to be the driven dipole voltage and V 2 zero, we 
eas'ily obtain the self admittance of dipole 1 and the mutual 
admittance (since Y 12 = Y^ ) between the dipoles by dividing the 
appropriate complex Fourier transforms of V,, I 1f and I 2 . 

The Moment Method results were obtained using the Electromagnetic 
Surface Patch Version 4 [7] computer code. The wire radius for 
the Moment Method calculations was taken as 0.281 cm, providing 
the same cross section area as the 0.5 cm square FDTD cells. 

While the FDTD calculations should be valid up to approximately 3 
GHz based on having 10 FDTD cells per wavelength, the thin wire 
approximation for the Moment Method code becomes questionable at 
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approximately 1 GHz and this was taken as the upper frequency 
limit for comparison of results. 


Ill Results 

Figure 2 shows the Gaussian pulse voltage applied to the 1 cell 
gap at the center of the longer, driven dipole. Figures 3 and 4 
show the current flowing in the center cell of the driven and 
passive dipole respectively. All are plotted on the same time 
scale, corresponding to 8,192 time steps. This calculation 
required approximately 7 hours on a 25 MHz 486-based personal 
computer . 

Figures 5-7 show the magnitude of the Fourier transforms of the 
voltage and current results of Figures 2-4. The current results 
indicate the complicated frequency domain behavior of the coupled 
dipole system. 

The self admittance was obtained by dividing the complex Fourier 
transform of the driven dipole current by that of the Gaussian 
voltage pulse at each frequency. The results are shown in 
magnitude and phase in Figures 8 and 9 and compared with ESP4 
Moment Method results. Considering the differences in how the 
feed region is modeled (a 1 cm gap in the FDTD calculations vs an 
infinitesimal gap in ESP4) the agreement is quite good. 

The mutual admittance was obtained in a similar manner, dividing 
the complex passive dipole current by that of the Gaussian pulse. 
The results are shown in Figures 9 and 10. Again the agreement 
is quite good considering the different approximations and 
assumptions made in the FDTD approach relative to the ESP4 
computer code. 


IV Conclusions 

The capability of the FDTD method to predict mutual coupling 
between antennas was demonstrated. The test case was two 
parallel wire dipoles of different lengths, with one driven by a 
Gaussian pulse. The complex self and mutual admittance results 
obtained using FDTD showed good agreement with results obtained 
using the Method of Moments. 
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